Comparison study of phase-field and level-set method for three-phase systems including two minerals

We investigate reactive flow and transport in evolving porous media. Solute species that are transported within the fluid phase are taking part in mineral precipitation and dissolution reactions for two competing mineral phases. The evolution of the three phases is not known a-priori but depends on the concentration of the dissolved solute species. To model the coupled behavior, phase-field and level-set models are formulated. These formulations are compared in three increasingly challenging setups including significant mineral overgrowth. Simulation outcomes are examined with respect to mineral volumes and surface areas as well as derived effective quantities such as diffusion and permeability tensors. In doing so, we extend the results of current benchmarks for mineral dissolution/precipitation at the pore-scale to the multiphasic solid case. Both approaches are found to be able to simulate the evolution of the three-phase system, but the phase-field model is influenced by curvature-driven motion.

typically be caused by minerals dissolving or precipitating, in case the alteration of the mineral layer is non-negligible. When the evolution of the fluid-solid interface depends on a solute concentration transported in the fluid, the evolution is not known a-priori and we have a free-boundary problem. To capture such evolving pore-scale geometries, most commonly level-set or phase-field methods are applied. The level-set method captures interfaces separating different phases as lower dimensional submanifolds within the computational domain, while the phase-field approach utilizes smoothed indicator functions for phase separation.
A comparison of level-set and phase-field methods in the context of precipitation/dissolution of a solute is found in [34]. Likewise, various approaches including levelset and phase-field methods were recently investigated in benchmark scenarios for mineral dissolution in [20]. Furthermore, both geometry evolution methods were compared regarding their respective strengths and weaknesses for the simulation of multi-phase flow problems on fixed domains in [2]. In contrast to these considerations, we apply the two approaches and compare them in the situation of a threephase system. Extending a phase-field model to the ternary case is feasible through generalizing the usual doublewell potential to a triple-well. Applications to three-phase flow [4,5] and two-phase flow together with one evolving solid phase [25,26] have earlier been derived and analyzed, and we here extend these approaches to the case of one fluid phase and two competing minerals. In this setting, we compare the level-set and phase-field approaches by means of simulation scenarios, but also in terms of sharp-interface limits. For modeling issues related to the case of one fluid phase and two mineral phases we refer to [13] in the context of levelset and to [15] in the context of phase-field approaches.
Contrary to the two-phase system, in which only one fluid phase and one mineral phase are present, interfaces are likely to develop high local curvature in three-phase systems even for regular initial conditions. Thus, their numerical treatment poses additional challenges in particular in the proximity of possible triple points. For the level-set approach, we apply the Voronoi Implicit Interface Method (VIIM) as presented in [10,13,29]. This method has been successfully applied to several physical problems such as curvature flow, multi-phase fluid flow and foam dynamics [29]. The phase-field equations do not need any particular tracking of interfaces and triple points, but the Allen-Cahn phase fields, which are considered here, include implicitly curvature-driven motion of the diffuse interfaces [3], which can lead to unwanted effects [30].
Inspired by [9,20], three scenarios with increasing complexity are studied. First, the ordinary differential equations (ODEs) related to a dissolution-precipitation reaction system are examined together with the corresponding evolving pore-scale geometry. Thereafter, the model is extended to additionally cover diffusive transport processes. Finally, single-phase fluid flow and advective transport are included into the model.
For all three scenarios, characteristic quantities such as the volume occupied by the individual mineral phases, the mass of the mobile species and the mineral's surface area are investigated and compared. Moreover, the time-dependence of the corresponding effective diffusion and permeability tensors is additionally evaluated.
The paper is outlined as follows: In Section 2, we introduce the chemical reaction system under investigation as well as the transport and flow model in their level-set and phase-field formulation. Additionally, the sharp-interface limits are discussed briefly. This is followed by a description of the numerical methods used for both approaches in Section 3. In Section 4, we compare and discuss the simulation results of all three scenarios including characteristic quantities and the evaluation of the effective diffusion and permeability tensors. Finally, in Section 5, we point out directions for further research.

Mathematical modeling
For an overview of quantities used in this paper, see Table 1. We consider a three-phase system in the two-dimensional, level-set function q advective velocity LT −1 ρ D , ρ P mineral densities NL −3 σ specific surface area L −1 V volume L 3 v n normal interface velocity LT −1 W well potential ξ diffuse-interface width L rectangular domain Ω = (0, 1)×(0, 1) with exterior boundary ∂Ω. More precisely, the domain Ω is time-dependently decomposed into the fluid domain Ω f (t) and the solid part Ω s (t) for all times t in the time interval (0, T ). At initial time t = 0, the solid Ω s (0) comprises two different mineral phases denoted D and P, cf. Figure 1. The fluid-solid interface (interior boundary of Ω f (t)) is disjointed into Γ int,D (t) and Γ int,P (t) accordingly. Furthermore, we denote the interface separating the two minerals by Γ int,M (t).
The fluid-solid interfaces and mineral phases individually evolve with time according to heterogeneous reactions which are specific to the two minerals, cf. Section 2.1 below. The mobile reaction partners may potentially be transported within the fluid domain Ω f (t) by molecular diffusion and advection.
In this paper, we investigate the following three situations of increasing complexity: First, the reaction system for precipitation and dissolution of the mineral phases is investigated. Disregarding diffusive and advective transport, the temporal evolution of the mobile species' concentrations within the fluid domain is described using ordinary differential equations (ODEs). Accordingly, the interface parts Γ int,D and Γ int,P each move with spatially uniform but time-dependent velocity. Since this approach does not resolve the potential spatial distribution of the concentration fields, it is only a valid approximation for regimes featuring strong diffusion and slow reactions. In the single-mineral case, such a simple setting with uniform interface velocity would not require a level-set nor a phase-field formulation, as an ODE for the mineral radius would be sufficient to describe the evolving geometry. Such a simplified approach was for instance first investigated in [21]. However, due to the interactions between the two competing minerals, level sets or phase fields are needed to resolve the geometry in the more complex two-mineral setting.
Secondly, we additionally resolve the spatial distribution of the concentrations. To this end, we consider a transport equation for the mobile chemical species including molecular diffusion and describe the chemical reactions as surface reactions on the two distinct fluid-solid interfaces (level-set approach) or highly localized volume reactions neighboring phase boundaries (phase-field approach). As the resulting chemical reactions are no longer uniform with respect to space, the geometry evolution becomes more challenging and results in the formation of complex mineral shapes. This again requires the usage of sophisticated methods in terms of modeling and numerics, also for the single-mineral case, as provided by means of a level-set or phase-field approach [21,22].
In our final simulation scenario, we additionally include the transport of the mobile species by advection. The related division of the interfaces into upwind and downwind parts with respect to the flow direction further increases the complexity of the solid-solute interaction and the evolution of the mineral shapes. Along these lines, we underline the capability of the level-set and phase-field methods.
Finally, the time dependence of the effective diffusion and permeability tensors is evaluated for all three scenarios. Based on upscaling theory, auxiliary cell problems are solved and their (flux) solutions are averaged as described in Section 4.2.2 for the level-set and phase-field approach.
In all settings, the chemical reactions drive the evolution of the pore-space geometry into an equilibrium state. In order to capture this evolution, level-set and phase-field methods are applied and compared. For a detailed description of the modeling and implementation details see Sections 2.2 and 3.1 for the level-set approach and Sections 2.3 and 3.2 for the phase-field approach, respectively.

Chemical reactions
In this section, we present the chemical setup of concern. We consider a reaction system describing dissolution/precipitation, e.g. the carbonation of silicates, as introduced in [9]. This involves the dissolving mineral D and the precipitating mineral P as well as mobile species A, B, C with reaction paths The reaction kinetics are chosen according to the classical law of mass action for a one-sided chemical reaction by including the back reaction via an equilibrium condition, cf.
(2), [9]. Since the chemical system in our simulations will be deflected in such a manner that both reaction paths in (1) effectively proceed in one direction from left to right, it is meaningful to refer to D as the dissolving mineral which will either dissolve or reach an equilibrium state when the net reaction rate is zero, while P denotes the precipitating mineral which will either precipitate or encounter zero net reaction rate. As this paper is mainly focused on the geometry evolution, we refrain from incorporating textbook data for a specific chemical system and formulate the equations without units. Nevertheless, dimensions are provided in order to clarify the physical meaning of the presented quantities, such as length L, mass M, number of particles N and time T. We use the following volumetric reaction rates [LT −1 ] for the concentrations of the mobile with reaction constants k D , k P [LT −1 ] and equilibrium constants K D [-], K P [L 6 N −2 ]. Throughout this paper, precipitation of a mineral is assumed to occur on the surface of that mineral only. Hence, we do not consider nucleation. Within the following sections, we introduce the three different modelling scenarios of increasing complexity each presented using a level-set and phase-field formulation.

Level-set formulation
In the level-set framework, we represent the geometry contained within the domain Ω using a real-valued function Φ : R 2 → R, the level-set function [23,29,32]. For further details we refer to Section 3. Accordingly, the time evolution induced by a space and time dependent normal velocity v n is implemented by solving the following advection equation for the level-set func- where the parameter v n : R 2 → R, [LT −1 ] prescribes the normal velocity with which level sets are transported. For our application, we define: Note that the normal velocities are defined only at the interfaces themselves and require to be suitably extended to the whole domain Ω, see Section 3.1 for details. Level-set models as derived here are well suited for formal upscaling into a two-scale model using periodic homogenization [13,21].

ODE model
In the ODE case, the chemical reactions presented in Section 2.1 are modeled by a system of three coupled ODEs, one for each solute concentration. This approximation is suitable if negligible spatial variation in the concentrations is expected.
In order to obtain molar reaction rates from the volumetric quantities f D , f P [L 3 (L 2 T) −1 ], we multiply by the molar density ρ D , ρ P [NL −3 ] of the reacting minerals. Assuming a closed system, concentrations additionally evolve driven by displacement as a secondary effect. That is, precipitation leads to a shrinking fluid domain Ω f and therefore an increasing concentration of all solute species while dissolution decreases concentrations. Summarizing, the total change of the species' concentrations is then given by the two volumetric reaction rates on the different mineral interfaces multiplied by the difference in particle density between solid and fluid, each scaled by the respective specific reactive surface area for i ∈ {D, P}. Note that this quantity is 'specific' with respect to the time-dependent fluid volume. These considerations lead to the following ODE system for the concentrations of the mobile species A,B,C: The prescription of normal velocities according to (4) ensures the conservation of mass within our model. Furthermore, the level-set framework allows to conveniently obtain a piecewiselinear approximation of the phase-separating interfaces. As such, characteristic geometrical quantities such as phase volumes or interface lengths σ i (t) are simple to derive from the level-set function Φ at any time t, cf. Section 3.1.

Diffusion model
In contrast to the ODE case, we now include diffusive transport into our model, i.e. the vector of concentrations c comprising the three mobile species A, B, C solves the following diffusion equation within the fluid domain Ω f (t): using a scalar and uniform molecular diffusion D m > 0, [L 2 T −1 ]. For convenience only, we suppose that the same diffusion coefficient applies for all mobile species. In order to model the insertion or extraction of solute particles to or from the fluid domain due to the heterogeneous chemical reactions described in Section 2.1, we impose the following flux conditions at the two distinct parts Γ int,D (t), Γ int,P (t) of the interior boundary: x ∈ Γ int,D (t), x ∈ Γ int,P (t), using the normal interface velocities v n,i , i ∈ {D,P} from (4). As such, the solute species concentrations are coupled to the geometry via time-dependent interior boundaries. The back coupling from the concentrations to the geometry is again realized by equations (4).

Flow model
In order to additionally account for an advective flux, we consider in Ω f (t) with advective velocity field q [LT −1 ] in our final example. For a matrix in column representation A = a (1) , . . . , a (n) we define ∇ · A = ∇ · a (i) i as the column-wise divergence. The boundary conditions supplementing (9) are again (7)- (8) as no advective flux is considered traversing the interior boundary, cf. (11). The velocity field is given by the solution of the stationary, incompressible Stokes problem with viscosity μ > 0, [ML −1 T −1 ] and the related pressure field p, [ML −1 T −2 ]. In this example, boundary conditions are chosen as follows: using Γ int (t) = Γ int,D (t) ∪ Γ int,P (t). Note that using the above no-slip condition at the solid-fluid interface, potential non-zero fluid velocity in the normal direction at the boundary is disregarded. Depending on the density differences between solute and minerals, some small nonzero velocity could appear in the normal direction [21]. Neglecting this non-zero normal component is justified in most applications since the solid-fluid interface velocity is small compared to the fluid flow as discussed in [18].

Phase-field formulation
The phase-field formulation does not consider sharp interfaces between the different physical phases but instead models them using regularised characteristic functions. These so-called phase-field variables φ i [-] (for i = 1, 2, 3) are evolved such that they always sum to 1. With these phase fields, we aim to approximate the equations of the level-set formulation in order to describe the same physical system as presented in Section 2.2. In our setup, φ 1 corresponds to the fluid phase, φ 2 to the mineral phase D, and φ 3 to the mineral phase P. Within each bulk phase the corresponding phase field is equal to (or close to) 1, with all other phase-field functions being equal to (or close to) 0. Near the interface of two phases the respective phase fields transition smoothly between 0 and 1. The width of this diffuse-interface regions is regulated by the phase-field parameter ξ > 0 [L] and a limit process for ξ → 0 in Section 2.3.4 is able to recover a sharp-interface formulation corresponding to the one used in the level-set model above but with an additional term corresponding to curvaturedriven interface evolution. The position of the interface can be approximated by the contour line of φ i = 0.5 or be recovered using a Voronoi method as for the level-set formulation. However, the simulation itself does not require such a reconstruction. Instead, the diffuse-interface region is captured by indicator functions γ ij and in these regions the chemical reactions impact the evolution of the phase-field variables. The variables for the dissolved concentrations are defined on the entire domain Ω and the boundary flux conditions (7), (8) are integrated into the conservation equations.
Thus the interfaces need not be tracked and no special consideration is necessary for the description of contact lines and triple points. The phase-field model is well suited for upscaling into a two-scale model using periodic homogenization [6,25].
We here use a phase-field model that is an extended version of the two models presented in [6,25]. In [25] a ternary phase-field model (two fluids and one mineral) is considered, and our phase-field evolution equations resemble the model there. The model in [6] considers only two phases (one fluid and one mineral), but also allows flow in the fluid phase. Hence, our extension to flow in Section 2.3.3 builds on this model.
The phase-field variables are evolved using the Allen-Cahn formulation [3]. In addition to a diffusive term this includes derivatives of a fourth order polynomial called the triple-well potential: Derivatives of this potential enforce the phase fields towards the values 0 and 1, while the diffusive term ensures a smooth, but steep transition between these values. The interface locations are captured using the indicator functions which attain their maximum for φ i = φ j = 0.5. These introduce a dependence of the phase fields on the chemical reactions by adding source termŝ Here, f D and f P are the respective reaction rates (2). With diffuse interface width ξ and diffusivity parameter ω [L 2 T −1 ] the evolution equations are (i = 1, 2, 3)

ODE model
As for the sharp-interface description with the level-set framework, the dissolved concentrations and thus the volumetric reaction rates and their contributions to the evolution of the geometry are first assumed to be independent of the spatial variable x. These concentrations could be evolved using the same system of ODEs (5), altering only the method of determining the specific surface area. In the phase-field model the interfaces are not tracked directly so the interfacial areas and the fluid volume would be computed as However, in the phase-field model the shift of the interfaces is not solely driven by the chemical reaction, but also by curvature of the interface. This would only impact the above system of ODEs by changing the fluid volume, altering the change to dissolved concentrations slightly. The loss of mineral volume and the corresponding change in dissolved species would not be directly captured in this case. Instead of additionally approximating the geometric changes due to this curvature effect, we prescribe conservation of species in the entire domain rather than only within the fluid volume. Hence, our system of ODEs for the evolution of the solute concentrations is ∂ ∂t where capture the stoichiometric coefficients of the chemical reactions and molar densities ρ D , ρ P of the respective minerals.
When ∂ ∂t V D = σ D V f f D and ∂ ∂t V P = σ P V f f P , this setup is equivalent to the system of ODEs used in the level-set formulation.

Diffusion model
With spatially resolved concentrations transported by diffusion, we adapt (6) to account for the phase distribution. Since the phase-field model is defined in the entire domain, we add terms accounting for the species bound in the minerals instead of using a boundary flux condition for the fluid-solid interfaces. The equations are given as

Flow model
To include fluid flow and advection of dissolved species in the phase-field model the (Navier-)Stokes (10) are modified to account for the distribution of phases, similar as done in [6]. The existing terms receive factors of the phasefield φ 1 while two new terms are added to account for moving interfaces and ensure vanishing velocity q inside the mineral phases: The function g(φ 1 , ξ) enforces q = 0 in the sharp-interface limit and is chosen as The main properties of this function are g(1, ξ) = 0 and g(0, ξ) = Kξ −1 > 0, such that for φ 1 = 1 the stationary Stokes (10) is obtained, while for φ 1 = 0 only q = 0 remains. As the term g(φ 1 , ξ)q serves to enforce a disappearing velocity inside the solid phases, the parameter K controls how strongly this is enforced in the phasefield simulation, since φ 1 can be slightly positive inside the minerals.
In addition to diffusive transport the dissolved species c(t, x) are now advected by the velocity q:

Sharp-interface limit
For the limit of ξ → 0 a sharp-interface model can be recovered from the phase-field formulation. The steps recovering the sharp-interface limit build upon the ideas in [8] and are analogous to the limits derived for the corresponding models in [6,25], and are therefore only shown briefly.
The Allen-Cahn model used in this comparison introduces an additional curvature effect to the interface velocities, but aside from that the corresponding sharp-interface model is obtained. To show this, the main idea is using matched asymptotic expansions; outer expansions are used to recover the bulk phases, and inner expansions for behavior near the diffuse interface. These are then matched through matching conditions [8].
Assuming outer asymptotic expansions inserting them into the phase-field evolution (15) and gathering coefficients of powers of ξ , one obtains the leading order equation The stable minimizers of the potential W (φ) are the unit vectors e 1 , e 2 , e 3 ∈ R 3 , corresponding to the three bulk phases.
Inside the fluid domain the constant value φ 1 = 1 is obtained. Considering this value for the phase field and inserting expansions for φ out , c out and q out into the flow and transport equations (20), (19) and (23), the limit process of letting ξ → 0 recovers the sharp-interface versions (10), (6) and (9), respectively. Inside the solid phases φ 1 = 0 and the equations simplify to 0 = g(0, ξ)q, recovering the desired no-flow condition.
The boundary conditions are recovered from the same equations using so-called inner expansions. Between the bulk phases designated by the stable solutions φ ∈ {e 1 , e 2 , e 3 } there are transition zones where two or all three phase-field functions are positive. At such transitions between two phases we can define an interface Γ d as the 1/2-level-set of the involved phase fields. Along this interface we define our unknowns in curvilinear coordinates (r, s), where s describes the position along the interface, and r indicates the signed distance from the interface in the direction of the outer normal ν ξ , x = y ξ (t, s) + rν ξ (t, s). Introducing the scaled variable z = r/ξ , we consider the asymptotic expansions in ξ : with corresponding expansions for c and q. These are combined with rewritten derivatives and matching conditions. For outer expansions such as φ out 0 and fixed t, s, z the limit for ξ → 0 for z > 0 respectively z < 0 are written as φ out 0 (t, y 1/2± ). The matching conditions relate limits of inner expansions for z → ±∞ to the outer expansions at the interface, e.g.
Together this allows recovery of the interface conditions [8].
The leading order terms of the phase-field equations at the transition zone between φ 1 and φ 2 , which corresponds to the interface Γ int,D (t), are where v n,0 is the local normal velocity of the interface, κ 0 the curvature of the interface, and From (28) one arrives at an equation describing the shape of φ in 1,0 , while (29) yields the boundary condition v n, Analogously, one obtains the boundary conditions for the other two interior interfaces where Γ int,M (t) denotes the interface separating both mineral phases, cf. Figure 1. These boundary conditions differ from the sharp-interface velocities given in (4) and introduce a curvature-driven movement. The strength of this effect is limited by the diffusivity parameter ω but it cannot be chosen arbitrarily small. This parameter controls how quickly a reasonably sharp interface is enforced, hence a small value of ω allows the reactions to cause an overly diffusive transition zone.
Remark 1 Note that in the recent publication [1], a different scaling of the Allen-Cahn equation is suggested to avoid curvature-driven motion in the sharp-interface limit. Instead of (15), a form corresponding to is analyzed. As shown in [1], the curvature-driven motion of the interface disappears at the dominating order of ξ . However, the analysis in [1] is performed for a case without chemical reactions. When the reaction ratesf i are included, the outer expansions can no longer recover stable solutions corresponding to the three phases. Hence, such an approach would either require no chemical reactions, or a further reformulation of the reaction rates such that the sharpinterface limit can be recovered.
Remark 2 For a system with two phases, [35] presents a modified Allen-Cahn equation, where an extra term is added so that the curvature contribution in the limit ξ → 0 is cancelled. A naive application to this three-phase problem has the desired effect at the interfaces between two phases, but also affects the evolution of the triple point. Without an analytical consideration such as done for the unmodified equations in [7], the precise effect on the limit behavior is unclear. However, for the phase-field equations with ξ > 0, the introduced term violates the conservation of φ 1 + φ 2 + φ 3 = 1 near triple points.
Inserting the inner expansions into the remaining model equations yields the desired boundary conditions (11) for flow as well as (7) and (8) for transport. The leading order term of the continuity equation, is integrated over R with respect to z and applying the matching condition (27) yields, at the fluid-solid interfaces, where y 1/2− corresponds to a point x on the interface (φ 1 (x) = 0.5). Here ν 0 denotes the first order term of the interface normal. From the leading order term of the momentum equation, one can then obtain the desired boundary condition q out 0 = 0. Evaluating the transport equations near the interface yields the following equations. The leading term is After integration with respect to z and using φ in 1,0 > 0 as well as the matching condition (27) one arrives at namely that c in 0 does not depend on z. Considering the next order terms O(ξ −1 ), along the interface Γ int,P (t) where φ 2 = 0 and applying (32) and (36) yields the equation Finally, integration with respect to z and application of the matching conditions yields the boundary condition on the interface Γ int, Analogously, the corresponding condition with b D instead of b P is derived at the interface Γ int,D . While this recovers the sharp-interface conditions (7) and (8), the interface velocity v n differs between the two models in accordance to (31).

Discussion of level-set and phase-field formulation
In contrast to the phase-field model, in which the chemical reactions occur as right hand sides, cf. Section 2.3, the chemical reaction enter the level-set model as boundary conditions. This is due to the fact that the boundary region of phase fields has a positive volume while the level-set interfaces are of codimension one. As a remark, note that using the level-set method, the specific surface areas σ i (t) must be constructed from the level-set function at each time-step (for more details see Section 3.1), whereas in the phase-field model this quantity is given by a simple integral of phase-field functions, cf. (16). Furthermore, the level-set approach requires a complex numerical scheme (see Section 3) due to the need of reconstructing the actual interfaces at every time-step. The phase-field equations, however, can be solved using standard schemes, and the phase-field variables can be directly incorporated into the equations describing chemistry and transport. We note that both modeling approaches require the choice of artificial parameters influencing the solution quality. For the level-set approach, only two adjustable parameters are present. > 0 is used in order to control the size of the area where a Vonoroi reconstruction of the interfaces is applied while the frequency of reinitialization poses the other degree of freedom. Details on the effect and practical choice of these parameters are given in Section 3.1. The phase-field model introduces several parameters to deal with the diffuseinterface behavior. As seen in Section 2.3.4, the expected sharp-interface model is captured as ξ → 0, except for an additional curvature-driven motion. The role of the parameters ξ, ω and K, which the phase-field model relies on, and the used numerical values are discussed in Section 3.2. For the level-set as well as the phase-field approach, parameters will be chosen specifically to minimize approximation errors in interface position and the influence of curvature terms on a given mesh.

Numerical methods
In this section we provide detailed information on the numerical methods used for geometry evolution, fluid flow and solute transport. A list of spatial discretizations and timeintegration schemes used for the three different simulation cases is provided in Table 2 summarizing the different numerical approaches for level-set and phase-field models. In order to ensure comparability of the simulation results between both geometry-capturing methods, solution algorithms are chosen carefully in conjunction with robust and well-established low-order discretizations.

Level-set implementation
In order to capture the geometry using sharp interfaces, we make use of generalized level-set methods. More precisely, the Voronoi Implicit Interface Method [29] (VIIM) is applied as three interacting phases need to be treated. As such, the interfaces are encoded in the level sets of a real valued function Φ. In conjunction with an indicator function χ : R 2 → {0, 1, 2} the total of three phases (0: fluid, 1: mineral D, 2: mineral P) are distinguished in our setting [29,32]. While the indicators serve the purpose of identifying the bulk of the different phases, the interfaces separating them are captured within level sets of Φ.
The Voronoi Implicit Interface Method constitutes a generalization of classical level-set methods in the sense that it does not rely on the neighboring phases to be distinguishable by the sign of the level-set function. Clearly, this restriction would not allow for the formation of triple points.
Instead, VIIM uses an -shifted version of an unsigned distance function d V with respect to the interfaces to encode their position. Accordingly, at the initial time we set Φ = − d V for a small positive parameter 1. As proposed in [29], we use the doubled mesh size = 2h. This choice corresponds to the smallest possible to ensure finite-difference stencils at points neighboring the 0-level sets of Φ to completely stay within a single phase and not cross interfaces. Therefore, a stable evolution via the level-set equation (3) including a correct assignment of each computational node to its phase is ensured. Also note that the 0-level set of Φ corresponds to the -level set of d V . As Φ has the signed distance function property in a neighbourhood of its 0-level set, application of equation (3) transports the 0-level set in a numerically stable manner. In order to recover the position of the actual interfaces (corresponding to d V = 0 ), a Voronoi reconstruction step is performed. By having chosen the minimal reasonable , approximation errors within the reconstruction procedure are minimized.
We note that VIIM, like many other level-set methods, is unable to precisely conserve the mass of each single phase as we will discuss in more detail in Section 4. In the recent literature, several approaches are available to counteract this phenomenon beyond application of higher-order discretizations or mesh-refinement. In [28], a volume-reinitialization scheme is introduced actively correcting the mass error introduced by the level-set evolution. Alternatively, adjusted level-set equations are available using normal interface velocity corrections to improve mass conservation, cf. [33]. However, such approaches are typically developed to simulate incompressible multi-phase fluid flow where additional regularizing curvature-terms are present. These terms are not only nonphysical in our specific application but also require at least second-order spatial discretizations to evaluate properly. Since errors in mass are found to be reasonable low throughout our simulations in Section 4, we use plain VIIM as presented in [29] as a robust first-order accurate scheme. Despite the large number of available discretization options such as finite element or finite volume approaches, we further adhere to finite differences schemes as employed in [29] due to their straight-forward implementation on regular grids.
In our application, a movement of the interfaces in normal direction is induced by dissolution and precipitation reactions at the two different mineral-fluid boundaries, cf. (4). For this research, the numerical algorithm presented in [10] is therefore supplemented with the derivation of normal velocity field from the chemical concentrations according to (4).
Denoting the connected components of the 0-level set of Φ by Γ i , i = 1, 2, 3 our geometry evolution algorithm consists of the following steps, cf. [10].
1. Initialize geometry via Φ, χ 2. Calculate signed distance functions d i = d i (x) to the interfaces Γ i by solving the Eikonal equation using the Fast Marching Method [31].
5. Initialize velocity values v V n (γ ) in a neighborhood of the interfaces, calculated from concentration fields according to (4). 6. Calculate velocity extension v n by solving Evolve Φ for a small time-increment using the levelset equation (3) applying a suitable upwind scheme, e.g. [31].
Steps 2 through 7 are iterated until the final simulation time is reached. In order to improve stability and accuracy of the algorithm, periodic reinitializations of the level-set function by setting Φ = − d V are applied [29]. In particular, the implementation of step 5 poses several difficulties. For the initialization of the normal velocity field, nodes on both sides of the interface must be labelled according to the local concentrations. Furthermore, a change in sign is needed when crossing the interface as gradients of d V switch orientation. In order to meet these requirements the following strategy is applied: First, information from the level-set function and the indicators is used to identify all nodes neighboring a node which belongs to a different phase. This constitutes the set of all nodes that require an initialization value. Each node is then given a label from 1, . . . , 6 according to the following scheme: Initially, we assign each element to one of the three phase separating interfaces Γ int,M , Γ int,D , Γ int,P . Then, we further discriminate with respect to the side of interfaces the points are located.
Second, we identify all finite elements belonging to the fluid phase which at least feature one solid edge. We will further call them boundary elements. These contain the concentration values from which normal interface velocities are calculated. Matching initialization nodes and boundary elements by distance minimization, velocity data are prescribed according to the labels given before.
The implementation is performed within the MATLAB [19] compatible framework of RTSPHEM [11]. Besides finite differences Eikonal-and level-set solvers using the well-known first-order upwinding scheme by Rouy and Tourin [27], it includes mixed finite element solvers for transport equation (6) on irregular triangular meshes using a mass-conserving lowest-order Raviart-Thomas RT 0 /P 0 discretization. An exponential upwinding scheme is implemented to stabilize advection dominated transport. Additionally, our code allows for adaptive alignment of boundary elements' edges to a piecewise linear approximation of the interior boundary. More precisely, we track intersections of the interfaces with the edges of a fixed underlying mesh. Adding the intersection points to the set of nodes, an aligned mesh is generated by applying a Delaunay triangulation algorithm. Rewriting the inhomogeneous flux conditions (7) equivalently as source terms on the boundary triangles, the nonlinearities are resolved using Newton's method, cf. [12]. As the arising source term is highly localized, we iterate until the residual is decreased by at least six orders of magnitude. Moreover, time-stepping for the reaction PDEs is adaptively coupled to the CFL condition of the level set evolution.
In the advective PDE case, Stokes equations (10) are solved on the well-established and stable lowest-order Taylor-Hood elements P 2 /P 1 . In order to cope with the inherent saddle-point structure of the problem, we employ the iterative Uzawa algorithm [17]. Choosing a low relative tolerance of 10 −7 as the stopping criterion, high solution accuracy is ensured. Within each time-step of the overall solution algorithm, the Stokes flow field q is computed on the current geometry and subsequently passed to the transport equation (9) as a coefficient.

Phase-field implementation
The phase-field model is implemented in DuMu x [16] using a cell-centered finite volume discretization with two-point flux approximation for both the phase-field equations as well as the transport equations for the dissolved chemical species. In the flow scenario also the continuity equation is discretized using these control volumes and the pressure degrees of freedom are placed at their centers. Without fluid flow, the model equations are not very complex, with simple storage, flux and source terms. The system is solved using an implicit Newton solver with adaptive time-stepping according to the number of needed Newton iterations as indicator. The Jacobians are assembled using numerical differentiation.
To include tightly coupled fluid flow and advection, additional care should be taken. We solve for the velocity components on staggered grids with control volumes and degrees of freedom shifted by half a cell in the respective spatial coordinates. The degrees of freedom are placed at the centers of these new control volumes, which corresponds to the normal velocities at the centers of the faces of the original cells. In DuMu x this is implemented using multiple domains which are linked via a coupling manager [16]. The coupling manager shares the data between the socalled momentum and mass sub-problems and approximates values where there is no adequate degree of freedom available. The two sub-problems are not solved individually but assembled into a single matrix. With the fluid flow depending on the phase-field variables and their derivatives, the existing manager is extended to expose and approximate values at the desired points of the mesh.
As mentioned in Section 2.4, the phase-field model relies on several parameters that affect its numerical behavior. The parameter K controls how strongly a zero-velocity is enforced inside the solid phases and how the velocity develops within the transition zone. Investigating the flow inside the minerals and the velocity profiles near φ 1 = 1/2 for expected velocities can give a reasonable choice for this parameter. For the presented simulations in Section 4.5, K was chosen K = 10000, and n = 10. While a high value of K is required to prevent significant nonphysical velocities inside the solid phases, this term also suppresses flow in the diffuse interface, which can lead to an underestimation of permeability [6]. The choice of K should scale with the fluid velocities expected near the minerals. In Section 4.5 relatively high inlet velocities are used, which leads to quite high velocities along the mineral grain. When calculating the permeability in the cell problems in Section 4.2.2, much lower velocities are used and thus a value of K = 200 is applied here.
The phase-field parameters ξ and ω control the profile of the phase-field functions and the shape of the bulk phases. The diffuse-interface width ξ affects the steepness of φ i in the transition zones. Meanwhile ω balances the diffusive and potential-driven effects on the phase-fields against reaction and storage terms in (15), and also controls the impact of the curvature effect. A small value of ξ enables a better approximation of the interfaces but the choice of ξ is limited by the spatial resolution. This means, ξ should be large enough for the transition zone to be spread over multiple degrees of freedom and we use a value of five times the mesh size to resolve the interface. As such, ξ plays a similar role as in VIIM, cf. Section 3.1, and its choice is likewise constrained from below in terms of multiples of the discretization lengths.
A smaller value of ω allows for a temporarily diffuse transition zone and reduces the impact of interface curvature, which can lead to more overgrowth of solids by small mineral tendrils. This makes it a central parameter affecting the simulation quality. However, if the changes caused by the chemical reaction dominate the contribution of the triple-well potential in the phase-field evolution (15), the resulting variables φ are prone to attaining values in between 0 and 1 in a larger transition zone and no longer exhibit an interface character. These transition zones, however, stray significantly from a sharp-interface description and can cause additional challenges for the numerical simulation. In the presented simulations the diffusivity parameter ω is chosen as 2.5 · 10 −3 . Varying this value allows for finding a sufficiently small choice, which reduces curvature effects without losing cohesion of the diffuse interface. A small value of ω furthermore requires a sufficiently small value of ξ and thus a fine spatial resolution and increased computational effort.
In the equations for the conservation of dissolved species (19) and (23) as well as the modified incompressible Stokes (20) φ 1 enters as a multiplicative factor. To avoid degeneration of these equations, we instead use φ δ = φ 1 +δ, adding a small regularization parameter [25]. The exception being the function g(φ 1 , ξ) which is unmodified. This regularization slightly disrupts the conservation equations and the momentum (20) no longer collapses fully to 0 = q. In the presented simulations this value was chosen as 10 −10 for the ODE and PDE formulations and as 10 −6 for the flow model.

Simulations
In this section, we specify all physical parameters, initial and boundary conditions used in the simulation scenarios. We moreover define characteristic quantities such as the volume of the individual mineral phases, the mass of the mobile species and the mineral's surface area and recall how they are specified using the level-set and phasefield approach, cf. Section 4.1. Based on the characteristic quantities, the three simulation scenarios as presented in Section 2 are investigated. Moreover, based on upscaling theory, time-dependent effective diffusion and permeability tensors are additionally evaluated. We emphasize the similarities, but also outline the differences resulting from the two problem formulations by means of level-set and phase-field description. Note that the following simulations are performed in two spatial dimensions. In order to stress the physical meaning of all appearing quantities, we still refer to a mineral volume and a surface area.

Simulation setup
All calculations are performed using a regular 200x200 mesh. This corresponds to the finest resolution used in similar studies [20,34] and is therefore considered practically feasible in applications. Furthermore, grid convergence studies presented in the Appendix indicate a sufficient resolution and discretization order for our setup. Due to the different kinds of underlying discretizations used in this paper, the mesh resolution specified above only refers to the number of nodes which is the same in both our simulation approaches. The only exception is made in the discretization of the transport equation using the level-set approach and the framework provided by the RTSPHEM [11] toolbox, where nodes are adaptively added to the triangular mesh in order to align edges with the fluid-solid interface, cf. Section 3.1. To provide insights into the complexity of the different sub-problems arising in both level-set and phase-field approaches, we compare the number of degrees of freedom (DoFs) and computation steps in Table 3. Note that the solutions to all sub-problems obtained by iterative methods are computed to high accuracy (relative residuum smaller than 10 −6 or two subsequent residua with difference smaller than 10 −8 ) not to compromise the simulation results by numerical artefacts.
As it is apparent from Table 3, the level-set approach requires less DoFs for geometry description than the phasefield method due to the ability of encoding the whole information within a single function. However, flow and transport equations involve a larger number of Dofs in the level-set framework. This is essentially due to underlying triangular mesh (instead of quadrilaterals) to simplify mesh adaptivity, cf. Section 3.1. Moreover, we note that in the diffusive and flow scenarios the number of time-steps is lower in the phase-field simulations due to a fully implicit scheme simultaneously solving all sub-problems, allowing for larger time-steps. However, since the solute concentrations in the discretization of the ODE scenario are updated after each time-step, sufficiently small time-steps are required for both approaches. Finally, we note that the number of Newton-steps is significantly larger in the levelset approach in exchange for higher solution tolerances of the linear solver (less Krylov steps). However, accuracy is solely determined by the non-linear residuum.
At initial time t = 0, a circular solid inclusion is placed in the unit square Ω with midpoint (0.5, 0.5) and radius r = 0.2 for all three scenarios, see Fig. 1 for a to-scale visualization. For simplicity, we assume the minerals D, P to be arranged in two hemicycles of the circle. For the phasefield model this corresponds to initializing φ such that without contributions from reactions or curvature effects the initial conditions are close to stationarity. This is achieved by using the base kernel of  (20). Finally, we choose the following initial conditions for the PDE cases as well as for the ODE model Section 2.2.1 disregarding the spatial variable x. Note that for the phase-field model, these initial conditions are chosen for the entire domain Ω, but correspond (in the sharp-interface limit) to (41) as φ 1 c is the relevant quantity. According to (41), solute species B and C are in chemical equilibrium at the initial time. Due to the oversaturation with respect to solute A, mineral D will immediately start to dissolve according to reaction (1) and release B. The resulting oversaturation with respect to species B will then trigger the precipitation of mineral P. The final simulation times T for each individual simulation are chosen in such a way that the steady state is approximated to a good extent, i.e. no further qualitative change of the system is expected to occur beyond that time.
Finally, we specify details on the boundary conditions on the exterior boundary ∂Ω for the different scenarios. For the diffusive PDE system, we choose homogeneous Neumann boundary conditions at the exterior boundary ∂Ω which correspond to no-flux conditions: i ∈ {A,B,C} For the advective flow case, we implement the following conditions at the outer boundary for t ∈ (0, T ) corresponding to an inflow boundary on the left, no-slip at top/bottom and an outflow boundary at the right hand side of the domain. Note that the flow at the left boundary shows a parabolic profile as expected for a Stokes flow within a long pipe. This results in a total volume flux of In order to attain an equilibrium state over time, we adapt the boundary conditions of the transport equations at the inlet (42) in the following way: x ∈ Γ inlet , i ∈ {A,B,C} using the equilibrium concentrations calculated in Section 4.3.1. Thus, the inflow will flush the over-saturated fluid domain and push the system towards the equilibrium state. We note that the resulting equilibrium volumes differ significantly from the characteristics of the previous cases as the flux boundary conditions dynamically change the total mass of solute species contained in Ω.

Measures
In the following, we define the five characteristic measures to quantitatively evaluate and compare the quality of the performed simulations for the different approaches.

Direct measures
The first measure is the mineral volumes and their evolution of over time as regarded in [20]. In the level-set approach, this quantity is easily inferred from a linear interpolation of the unsigned distance function d V on the underlying grid, see Section 3.1. Using the phase-field method, the mineral volume is given as the integral of the respective phase-field function φ i , see Section 2.3. Given the setup presented in Section 4.1, the initial volumes of both minerals amount to 0.2 2 π 2 ≈ 0.0628. As a second measure, we compare the surface area of both minerals, cf. [20]. That is, we compute the length [L] of the interior interfaces Γ int,D , Γ int,P , cf. Figure 1, separating the fluid domain and the respective mineral. In the given context, the resulting quantity therefore equals the reactive surface area. In the level-set framework, a linear reconstruction of the interface is used to approximate its length. Suitable interface indicators are used to obtain the related quantity using phase fields, cf. (16). Given the setup presented in Section 4.1, the initial surface areas of both minerals amount to 0.2π ≈ 0.628.
The third measure is the conservation of mass with respect to each chemical species. Although the formulation of the reactive problem presented in Section 2 is analytically mass conservative, the numerical schemes may not be capable of preserving this property precisely. Furthermore, as indicated in Section 4.3.1, orbits in the system's phase space regarding different total masses M have a positive distance. Hence, the relative loss or gain in total mass is a useful indicator to assess the simulators' predictive power as already used in [13].
Taking the fluxes j i = D m ∇c i − qc i related to solute i, i ∈ {A,B,C}, at the inlet and outlet (Γ inlet , Γ outlet ) for the advective example into account, the total mass M B (t), M C (t) of species B,C in the PDE cases is given by (cf. (51)) In this expression, the first term accounts for the amount of species i being dissolved in the fluid whereas the second describes the amount ligated within the minerals. Last, the third term is related to mass exchange across the domain's boundaries. In the context of the phase-field model the first term is adjusted slightly, integrating instead over the entire domain and including a factor of φ 1 to account for the fluid phase.
Note that the related quantity M A (t) is not considered in this paper as the solid shares would need to be weighted with a negative sign compromising physical interpretability. As the fluxes are explicitly discretized in both numerical schemes (Section 3), total in-and outflow are simple to determine by integration in the level-set as well as phase-field framework. Volumes and integrated concentrations in (45) and (46) are derived similarly.

Effective measures
Finally, we consider two effective quantities derived from the geometrical setup. Assuming the domain of interest Ω to be a representative elementary volume of a larger scale porous medium, we can ask for the effective diffusion and permeability tensors of that respective medium. These quantities are of high importance concerning flow and transport properties on a macroscopic scale. Both tensors are derived solving different auxiliary PDEs (cell problems) using periodic boundary conditions on the exterior boundary. In the context of periodic homogenization, the following sharpinterface representation is derived, cf. [14] for the static case and [21] for time-evolving domains: The diffusion tensor is given as for i, j ∈ {1, 2} and Kronecker delta δ ij , where ζ j are the solutions of the elliptic problems for j ∈ {1, 2} and outer unit normal ν, denoting again the total interior boundary by Γ int (t) = Γ int,D (t) ∪ Γ int,P (t).
The permeability tensor is given as In case of the phase-field framework, the above equations change their form. To calculate effective diffusion and permeability tensors we also solve auxiliary cell problems incorporating the phase-field parameter, with periodic boundary conditions on the exterior boundary. However, in the phase-field formulation the domain is not split along an interior interface, hence the cell problems are solved not only in the fluid domain Ω f but in Ω. Boundary conditions on the interior solid boundary are hence not needed, as they are already incorporated in the phase-field formulation. The effective tensors are then calculated with the regularization factor φ δ [6] and are defined as, for i, j ∈ {1, 2} with ζ j and (ω j , π j ) solutions to [6] − ∇ · φ δ D(∇ζ j + e j ) = 0 inΩ, As the off-diagonal components of the effective tensors remained small and generally several orders of magnitude smaller than the diagonal components, we will restrict to reporting and discussing the evolution of the diagonal elements.

Comparison for ODE system
For the ODE case it is possible to theoretically deduce the system's long term behavior for several characteristic quantities. As such, we start by providing important analytical properties of the ODE system.

Theoretical considerations
In this section, we briefly discuss existence of solution to the ODE problem (5) introduced in Section 2.2.1 as well as stability and positivity of equilibrium points.
Assuming continuous dependence of fluid volume and interface lengths on time, system (5) admits a unique solution local in time by the Picard-Lindelöf theorem. In case those quantities and the concentrations are bounded from above and away from zero, the solution can be extended globally.
Next, we are concerned with the stability of equilibria. In the following, we approach stability by investigating the equilibrium points as a function of the system's invariants. Let us denote the volume of minerals D and P being present at time t by V D (t) and V P (t), respectively. Due to the conservation of mass, the following quantities are conserved over time, cf. (45), (46) in Section 4.2.1: Accordingly, we have the following function mapping G from the system's state space to a set of invariants and characteristics for an equilibrium state with M = (M A (0), M B (0), M C (0)): More precisely, G maps the current solute concentrations and mineral volumes to the total masses and interface reaction rates. Clearly, all equilibrium points for a given M are characterized by the preimage G −1 (M, 0, 0) assuming σ i (t) being bounded from below by a positive constant. As such, a possibly continuous G −1 would lead to a curve of equilibria points in the phase space, rendering linearized theory inconclusive due to a zero eigenvalue. Furthermore, all equilibria reachable from positive initial conditions are located within the positive octant as easily seen by application of the quasi-positivity theorem [24]. Evaluating the expression for the initial conditions given in (41) we find c A,eq = c B,eq = 1.4056, V P,eq = 0.1333, V D,eq = 0.0339, as an equilibrium state. Investigating the Jacobian ∇G at that point shows local bijectivity of G. Accordingly, the system is not asymptotically stable. This is expected as disturbances changing M cannot be compensated by the system due to conservation of mass. As such, discretization errors with respect to the geometry evolution will add up. This inherent property underlines the necessity for welldesigned numerical algorithms.

Comparison simulations
The simulation results for the ODE case as outlined in Sections 2.2.1, and 2.3.1 are depicted in Fig. 2 for different time-steps and in Fig. 3 the evolution of quantitative measures is shown until the final simulated time of T = 2. At final time of the level-set simulation, the concentrations of all solute species deviated by less than one per mille from their calculated equilibrium values (53). As such, the system reaches equilibrium to high precision. In the case of the phase-field model the solved system of equations corresponds to a modified sharp-interface formulation, with additional curvature-driven interface motion. The simulation approaches an equilibrium with constant curvature along each interface and dissolved concentrations slightly perturbed relative to the calculated equilibrium such that the curvature-driven motion and reactive effects cancel out. At final time of the phase-field simulation the concentration c C matches up to three per mille but c A is two percent lower and c B two percent higher. Figure 2 displays and compares the resulting geometrical configuration of the three-phase system for different simulations times and both approaches. More precisely, the minerals as obtained using the level-set method are shown in grey (P) and black (D). In the surrounding fluid domain, the (spatially independent, but time-dependent) concentration of solute species A is displayed. Finally, the phase-separating interfaces as obtained using the phase-field model are overlaid in white.
Overall, we observe a good match between level-set and phase-field simulation results. However, it is evident that the level-set shows overgrowth of the minerals, while the phase-field does not, see Fig. 2. Due to curvature effects, the phase-field model cannot resolve the corners of the dissolving mineral very well, and the initially straight interface between the two minerals becomes curved.
Although the shapes evolve slightly differently, a good match in characteristic measures as depicted for both levelset and phase-field solution in Fig. 3 is observed. In particular, the volumes of the two minerals almost perfectly coincide as seen in Fig. 3. Both predicted volumes of mineral P at final simulation time differ less than 0.3% from the analytical equilibrium volume, for mineral D the relative deviation is less than 2.5%. These disagreements are considered fairly small given that, starting from an initial volume of 0.0628, equilibrium volumes of 0.1333 and 0.0339, respectively, are targeted. As such, initial volumes are approximately doubled (mineral P) or halved (mineral D) in the course of the simulation. The evolution of the surface area is, except from the initial period, also comparable. Due to the large increase of mineral volume, the surface area of mineral P almost doubles within our simulation from 0.628 to 1.215 in the level-set simulation, see Fig. 3. Simultaneously, originating from the same initial value, the surface area corresponding to mineral D decreased by more than 35%. Approaching the equilibrium state, both modelling methods concur well with relative differences of 5% and 11% for minerals D and P, respectively.
The initially peaking surface areas for the phase-field simulation seen in Fig. 3 are caused by how they are determined from the phase-field variables φ. Without reactive contributions the phase-field variables develop a specific profile and the surface area is calculated as the integral of 4 ξ φ i φ j . During the early evolution in the ODE case, the reaction rates are high and the shape of the phasefield variables fail to keep the expected shape. The transition zones get drawn out and the changed profile across the interfaces causes the above integral to overestimate the interface length. As the phase-field ODE formulation does not depend on this explicit evaluation of the surface area, this does not have a strong effect on the further evolution of geometry and concentrations. We further note that the formulation used for the phasefield approach is perfectly mass conservative (up to 10 significant digits), while small deviations within 0.5% are seen for the level-set approach. As indicated by the grid convergence studies performed in Fig. 8, mass errors in the level-set framework decrease consistently with higher spatial resolution.
Finally, Fig. 3 shows an almost perfect match of the effective diffusion and permeability tensors for small times t ≤ 0.3. Across both directions and tensors, relative to that time according to the phase-field simulation. As the chemical reactions come to a standstill, the level-set data plateau whereas the phase field data further evolve due to curvature effects. For both effective quantities, the phase-field simulation reduces the distance between the measurements in x and y direction. Apparently, curvature effects diminish anisotropy over time.

Comparison for diffusion model
The simulation results of the diffusive PDE case as outlined in Sections 2.2.2, and 2.3.2 are depicted in Figs. 4 and 5 for both approaches. The concentration of solute species A is displayed according to the level-set simulation. Intermediate and a close to equilibrium state of both mineral phases are illustrated for the level-set method including also the phasefield solution as an overlay, cf. Figure 4. As both simulation agree on generating mineral overgrowth (in contrast to the ODE case), the modeling approaches recover important qualitative physical behavior of the underlying problem. This is a direct consequence of the increased and contrast rich interface velocities at the triple points. These again result from focusing reactive activity to the neighborhood of the triple points made possible by taking the spatial distribution of the solute species into account. Therefore, the difference in interface velocity between both fluid-solid interfaces is significantly higher, facilitating overgrowth behavior. As such, the system is naturally forced to develop and maintain higher interface curvatures than in the ODE case. However this increased curvature causes a stronger deviation between the two models.
As earlier, the phase-field formulation cannot properly resolve corners due to the curvature-driven movement of the interfaces. Therefore, mineral D is increasingly displaced by mineral P within the solid close to the interface. This behavior is also visible in Fig. 5, wherein the concentration fields of the mobile species as well as the area/volume of the mineral species are investigated with respect to time. Due to the no-flux exterior boundary conditions, the system approaches an equilibrium state for large times similar to the ODE case discussed in Section 4.3. In fact, both systems' equilibrium states are identical in terms of mineral volumes and solute concentrations due to the same total masses of species A, B and C. Yet, the rate of convergence is much slower due to the time consuming transport of solute species between the two different mineral interfaces, see Fig. 4. As the transport speed is governed by the concentration gradients, convergence to equilibrium is additionally decelerated. At the final time  Diffusion case: Evolution of mineral volume over time, reactive surface area and relative total mass calculated with level-set method (LS, dashed lines) and phase-field approach (PF, dotted lines). The last two pictures illustrate the evolution of diffusion and permeability tensors along the main axes over time T = 4, the volume average concentration of solute A is 1.56. Accordingly, the transition from initial chemical disturbances to the equilibrium state is already completed by 74%. Yet, within the final simulation time unit 3 ≤ t ≤ 4 (cf. Figure 4) an insignificant amount of geometry evolution is identified.
The overall mineral volumes in the two models evolve qualitatively similarly, yet slightly differently. More precisely, the volume of mineral P is predicted progressively higher by the phase-field than by the level-set method (19% at final time) and vice versa for mineral D (22% at final time). Given a total volume growth of mineral P of 72% (from 0.0628 to 0.1081) as predicted by the level-set simulation, these deviations are not negligible yet reasonably small. A similar conclusion is drawn for mineral D, which shrinks by 32% from 0.0628 to 0.0424 over simulation time. This effect is assumed to essentially result from the artificial displacement at the mineral-separating interface in the phase-field simulation. In conclusion, we observe significantly higher deviations between the two simulation approaches compared to the ODE case presented in Section 4.3.2, due to the more pronounced mineral overgrowth and the resulting higher interface curvatures.
The surface area evolutions predicted by both approaches are also quite similar. At final time, we observe a relative deviation in mineral P's surface area of 9%. This is comparable to the deviation of 11% measured in the ODE case. According to the level-set prediction, the mineral underwent a growth of surface area by 111% from 0.628 to 1.325, which is about 17 percent points more than in the ODE case. Due to the significantly decreased surface of mineral D by 69% (from 0.628 to 0.197), relative differences appear to be high. Yet, in absolute measure, they are comparable to the deviations measured for the precipitating mineral.
Again, the phase-field model conserves mass up to 10 significant digits, while the level-set model experiences an error of about 2%. As illustrated in Fig. 8, the error in mass conservation consistently reduces with increased resolution in the level-set simulations.
Since the precipitation of P is concentrated at the poles of the initially circular geometry, both modelling approaches agree that only a slight change of effective permeability and diffusivity in y-direction takes place over time, cf. Figures 5  and 4. The associate values at final time differ by 12% and 3% among both modeling approaches, respectively. Due to the resulting increase in vertical extend, both quantities decrease significantly with respect to the x direction by 44% (from 0.0329 to 0.0185) and 13% (from 0.778 to 0.677), respectively, according to the phase-field computations. At final time, the results obtained by level-set and phase field approach differ by less than 2% for both effective tensors in y-direction.  Evolution of mineral volume over time, reactive surface area and relative total mass calculated with level-set method (LS, dashed lines) and phase-field approach (PF, dotted lines). The last two pictures illustrate the evolution of diffusion and permeability tensors along the main axes over time

Comparison for flow model
For this final example, transport of solute species is subjected to an additional advective flow field according to models ontroduced in Sections 2.2.3, and 2.3.3. Figure 6 displays and compares the resulting geometrical configuration of the three-phase system for different simulations times and both approaches. As in the previous cases (Figs. 2 and 4), the minerals as obtained using the level-set method are shown in grey (P) and black (D).
In the surrounding fluid domain, the (spatially and timedependent) concentration of solute species A is displayed. Interfaces as obtained using the phase-field model are overlaid in white. Apparently, both the level-set and phasefield approach nicely agree on geometry evolution. The system has reached an equilibrium state at the final time T = 1.5 to a reasonable extent. Throughout the domain Ω f (1.5), a maximal deviation of 6.4% from the equilibrium values is measured across all solute species. Similar to our previous scenarios curvature effects are visible in the phase-field simulation close to the triple points. Due to the reduced amount of mineral overgrowth compared to Section 4.4, the implications are less severe. Figure 7 proves good agreement in mineral volume (relative deviations of 8% and 2% for minerals P and D) as well as surface area prediction (relative deviations of 11% and 10% for minerals P and D) among both simulation approaches. Still this variance is reasonably small compared to the loss of 22% and 27% in volume and surface area for mineral D and gain of 18% and 41% in volume and surface area for mineral P over simulation time according to the level-set approach. The phase-field method again achieved almost perfect mass conservation, while the level-set method's error is below 1.5%.
Due to the reduced difference between initial and final geometry in comparison to the diffusive PDE and ODE case, also the evolution of the effective tensors is less pronounced. As such, both permeability and diffusivity in y direction remain nearly constant over time. Here, relative deviations of 6% and 2% are observed at final time between both modeling approaches for diffusion and permeability, respectively. The evolution with respect to the x-direction appeared again to be more significant. A decrease of 11% (from 0.0329 to 0.0292) and 2% (from 0.778 to 0.759) are measured according to the phase-field simulation, respectively. At final time, both models agreed on the effective quantities with a deviation less than 5% with respect to the x direction.

Conclusion
As shown in Section 4, both the level-set model and phasefield model are able to simulate geometrical changes in a three-phase system involving two minerals and a fluid including the solutes taking part in heterogeneous reactions. Upscaled quantities like permeability and diffusivity predicted by the two approaches show a comparable evolution as the mineral shapes evolve. However, each of the two approaches inhibit strengths and weaknesses.
As it is apparent from the simulation results in Section 4, the phase-field model is capable of conserving mass up to 10 significant digits. Using a finite volume discretization the presented phase-field model can conserve mass almost perfectly, with minor losses due to the regularization of the conservation equations. In contrast, the level-set method used does not generally guarantee conservation of mass. Yet, with a maximal deviation of 2% throughout our experimental lineup, the loss/gain in mass in the level-set simulations is relatively low. Furthermore, grid convergence tests conducted in the appendix (Fig. 8) show a diminishing effect on higher resolution. This is due to the fact that the changes in mineral volumes V i are determined bẏ V i (t) = v n,i (t)|Γ int,i (t)| for i ∈ D,P using an explicit first-order time discretization. In theory, the error in mass is controlled linearly by the timestep size. Assuming piece-wise smooth interfaces, the error could be reduced in higher-order by applying a higher-order discretization scheme or additionally taking local curvature into account.
On the other hand, the results of Section 4 certify the level-set method to properly handle high curvature within the interfaces. Although the reconstruction procedure applied within VIIM (see Section 3.1) loses accuracy in these situations, the errors remain highly localized and do not propagate along the whole interface over time. As such, high interface curvatures along the static interface separating both solid minerals are recovered very well. In contrast, the phase-field model performs increasingly unsatisfactory close to the equilibrium as interface velocity becomes curvature dominated. In general, the curvature-driven motion of the phase-field model affects the simulation results. As remarked in Section 2.3.4, existing approaches for diminishing curvature-driven motion are not applicable for the current setup. Although one can choose the relevant parameters controlling the curvature-driven motions small, they cannot be chosen zero and are limited by the choice of the grid. As a fine grid is needed to resolve the diffuse interfaces, a natural extension would be using adaptive grid refinement for the diffuse transition zones since one can generally accept a much coarser grid away from the interfaces. Furthermore, the phase-field model requires well chosen parameters, in particular the phase-field diffusivity ω. As described in Section 3.2 the phase-field parameters strongly influence the simulation quality, but for ω there is no simple way to predict a good choice.
Overall, the two presented modeling approaches and their implementations highlight the difficulty with capturing evolving interfaces attaining high curvatures. As the numerical experiments do not indicate any major difficulties for the level-set approach except from some mass loss/gain, the phase-field approach is influenced by the curvature-driven motion. Further work is required to find a suitable strategy to eliminate curvature-driven motion and determine suitable parameters in the phase-field approach. Phase-field: Simulation results for different spatial resolution and choice of the phase-field parameter ξ in case of the diffusive case (top), equilibrium flow case (bottom). Volume and surface area of mineral P as well as mass conservation with respect to solute B are presented already mentioned, a smaller choice of parameters increases precision of the geometry evolution. Yet, a lower bound is given by some multiple of the discretization length. As such, we investigate both the impact of mesh refinement as well as parameter reduction , ξ on constant meshes. The results of our convergence studies are presented in Fig. 8 for the levelset model and in Fig. 9 for the phase-field approach. For clarity, the data at final simulation time T are additionally listed in Tables 4 and 5.  [9] to our attention. We further acknowledge the insightful discussions with Florian Frank and Helmut Abels.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data Availability
The codes used in this paper will be made available within their respective frameworks of DuMu x [16] and RTSPHEM [11] upon publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.