Upscaling of a Cahn–Hilliard Navier–Stokes model with precipitation and dissolution in a thin strip

Abstract

We consider a phase-field model for the incompressible flow of two immiscible fluids. This model extends widespread models for two fluid phases by including a third, solid phase, which can evolve due to e.g. precipitation and dissolution. We consider a simple, two-dimensional geometry of a thin strip, which can still be seen as the representation of a single pore throat in a porous medium. Under moderate assumptions on the Péclet number and the capillary number, we investigate the limit case when the ratio between the width and the length of the strip goes to zero. In this way, and employing transversal averaging, we derive an upscaled model. The result is a multi-scale model consisting of the upscaled equations for the total flux and the ion transport, while the phase-field equation has to be solved in cell problems at the pore scale to determine the position of interfaces. We also investigate the sharp-interface limit of the multi-scale model, in which the phase-field parameter approaches 0. The resulting sharp-interface model consists only of Darcy-scale equations, as the cell problems can be solved explicitly. Notably, we find asymptotic consistency, that is, the upscaling process and the sharp-interface limit commute. We use numerical results to investigate the validity of the upscaling when discontinuities are formed in the upscaled model.

1. Introduction

Multi-phase flow and reactive transport in porous media are encountered in many important fields, including geological ${\rm CO}_2$ sequestration, geothermal energy, groundwater management, oil recovery and ion exchange in fuel cells. While the modelling of multi-phase flow is itself a challenging task, the examples given before have in common that the solid matrix of the porous medium can change in time due to processes such as precipitation or dissolution, which, in turn influence the flow behaviour.

Another common point of the processes mentioned before is that they take place in a porous medium. In this case, two different length scales are encountered. The first is the pore scale, where each phase (solid, or fluid) can be identified clearly, occupying certain positions in well-defined volumes. The second is the so-called Darcy scale, which is used in most situations of practical relevance, and where averaged quantities are used to describe the behaviour of the system.

Such processes can be modelled at different scales. When formulated at the pore scale, the models are capable of describing the detailed processes accurately. On the other hand, they are defined in a highly complex domain, the union of the pores of in the porous medium, and this makes such models difficult to use for real-life applications. Instead, Darcy-scale models are formulated without taking into account the detailed behaviour of the system at the pore scale, employing constitutive relationships that are stated directly at this larger scale. Therefore, one may say that Darcy-scale models are suited for practical applications, but are missing the accuracy of the pore-scale models. In this context, upscaling is a natural way to derive mathematical models that, on one hand, can be used for practical applications, and, on the other hand, do incorporate accurately the processes taking place at the pore scale. We refer to Dentz et al. (2011) for an overview of reactive transport models in porous media.

In detail, we are interested here in the situation where two immiscible fluid phases occupy the pore space of a porous medium. One fluid phase contains ions that can precipitate at the fluid–solid interfaces. This leads to the formation of a precipitate layer at the pore walls, which reduces the space available for the fluid. The reverse process, that is the dissolution of the mineral phase into the fluid phase, is also allowed. In this case, the volume of the precipitate is reduced, while the volume available for flow is increased, and more ions are dissolved in the fluid phase. More precisely, the pore-scale model uses the conservation of mass, momentum and of the dissolved ions in each phase. The challenging aspect here is related to the fact that the spaces occupied by each of the two fluids, as well as by the mineral, can change over time in an a priori unknown manner. Therefore, the different phases are separated at the pore scale by free boundaries, which are unknowns in the model. After upscaling, these free boundaries translate into unknown Darcy-scale quantities such as fluid saturation, mineral concentration or the porosity and permeability of the medium. In particular, the latter become time- and space-dependent unknowns, satisfying evolution equations. This is not the case of commonly used Darcy-scale models, where, as mentioned, a given relationship between quantities at the Darcy scale is assumed (e.g. the Cozeny–Kármán relationship).

Different approaches have been proposed when developing mathematical models for applications involving free boundaries at the pore scale. For a simple geometry, which is basically a long, thin strip (in two spatial dimensions) or tube (in three dimensions) the free boundaries can be viewed as functions of one or two variables. In this sense we mention van Noorden (2009a) for a model describing precipitation and dissolution but for one fluid phase, which has been extended in Agosti et al. (2016), Bringedal et al. (2015), Kumar, van Noorden & Pop (2011) and Kumar, Wheeler & Wick (2013) to incorporate non-isothermal or mechanical effects, or in different flow and reaction kinetics regimes but still for the saturated, single-phase flow, and Mikelić & Paoli (2000); Mikelić (2009); Picchi & Battiato (2018); Sharmin, Bringedal & Pop (2020); Lunowa, Bringedal & Pop (2021) for unsaturated single-phase flow or two-phase flow models.

For more complex geometries, level sets can be employed to describe the evolution of the free boundaries. In this respect we refer to van Noorden (2009b), as well as to Bringedal et al. (2016); Schulz et al. (2017); Schulz (2019), all considering models for precipitation and dissolution in a water-saturated porous medium.

When applying any of both approaches mentioned before, one has to deal with (freely) moving interfaces. This makes not only the mathematical analysis, but also the development of efficient numerical scheme, a challenging task. Alternatively, one can use phase fields to approximate the interfaces between phases by diffuse transition zones with small positive width. The phase fields are smooth approximations of the indicator function of each phase. The evolution of the phase fields is usually derived as the gradient flow to a free energy and, in the limit case when passing the diffuse-interface parameter, one should recover the original, free boundary model.

Commonly used phase-field models involve either the Allen–Cahn equation (Allen & Cahn 1979) or the Cahn–Hilliard equation (Cahn & Hilliard 1958). While the Allen–Cahn equation is of second order and ensures that the phase-field indicators remain essentially bonded by zero and one, it is not conservative. Therefore, we focus here on the Cahn–Hilliard equation, which is of fourth order but conservative for the phase-field indicators.

Models coupling the Cahn–Hilliard equations and the incompressible Navier–Stokes equations have been developed for two fluid phases (Abels, Garcke & Grün 2012), three fluid phases (Boyer & Lapuerta 2006; Boyer et al. 2010) and more than three fluid phases (Boyer & Minjeaud 2014; Dunbar, Lam & Stinner 2019). For the description of fluid–solid interfaces, the Navier–Stokes equations can be solved in the fluid volume fraction and a velocity of zero is assigned to the solid phase (Beckermann et al. 1999; Sun & Beckermann 2004). Phase-field models are also used in Baňas & Mahato (2017), Bunoiu et al. (2020), Daly & Roose (2015), Metzger & Knabner (2021), Schmuck et al. (2012) and Schmuck et al. (2013) as pore-scale models for two-phase flows in porous media, and further Darcy-scale models are derived. Kinetic reactions at phase boundaries have been introduced in van Noorden & Eck (2011) and Redeker, Rohde & Sorin Pop (2016). The pore-scale model in Redeker et al. (2016) includes two immiscible fluid phases and a mineral one, but the fluid phases only move due to curvature effects. Also, the corresponding Darcy-scale model is derived by homogenization techniques. More recently, phase-field models that couple precipitation and dissolution with fluid flow have been developed in Bringedal, von Wolff & Pop (2020) (for one fluid phase, and for which the Darcy-scale model is derived), and Rohde & von Wolff (2021) for a two-phase flow.

The starting point in this work is the Cahn–Hilliard–Navier–Stokes model developed in Rohde & von Wolff (2021), which describes the processes at the pore scale. The aim is to derive an upscaled model corresponding to the Darcy scale. We consider the simplified geometry of a thin strip, and assume that the ratio of the width of the strip and its width is small. We employ asymptotic expansion methods that use this ratio as expansion parameter, and derive upscaled equations for transversally averaged quantities. In this respect, we follow the ideas in Bringedal et al. (2015), Kumar et al. (2011) and van Noorden (2009a) for one-phase flow including precipitation and dissolution effects at the pore walls, and Mikelić & Paoli (2000); Mikelić (2009); Sharmin et al. (2020) and Lunowa et al. (2021) for two-phase flow, all considering a thin strip or tube. Observe that the pore-scale models in these works mentioned above involve free boundaries. Instead, for the phase-field pore-scale model in Bringedal et al. (2020) describing the flow of one fluid phase but including precipitation and dissolution, a Darcy-scale model is also derived for a thin strip by transversal averaging, in comparison with the one obtained by homogenization in more general situations.

We recall that the geometry considered here is simplified, a long and thin pore. In this case, asymptotic expansion methods and transversal averaging are sufficient to derive the Darcy-scale models. For completeness, we mention that, for general geometries, different techniques may be needed. Restricted to processes involving free boundaries at the pore scale, and particularly to two-phase flows, with or without mineral precipitation or dissolution, both homogenization and volume averaging methods are suited for deriving Darcy-scale models. Homogenization is used in Bringedal et al. (2016), van Noorden (2009b), Schulz et al. (2017) and Schulz (2019) for models based on level sets, and in Bunoiu et al. (2020), Bringedal et al. (2020), Redeker et al. (2016) and Schmuck et al. (2013) for phase-field approaches. Alternatively, volume averaging methods are used e.g. in Bahar et al. (2016), Quintard & Whitaker (1988), Quintard & Whitaker (1994), Whitaker (1986) and Tartakovsky et al. (2007) to derive Darcy-scale models for two-phase flow or reactive transport in porous media. Finally, we also mention that Darcy-scale models for problems of the type discussed here can be obtained by the thermodynamically constrained averaging theory, as done in Gray & Miller (2005), Jackson et al. (2012) and Rybak, Gray & Miller (2015).

The main contributions here are threefold. First, starting from a pore-scale model, asymptotic expansion arguments are employed to derive a two-scale model for the two-phase flow in a porous medium, in which the dissolution and precipitation effects are taken explicitly into account. Compared with commonly used Darcy-scale models, instead of postulating relationships between the Darcy-scale quantities such as porosity and permeability, these are obtained here by solving (pore-scale) cell problems. Also, the situation considered here is more complex than in previous publications with a similar focus, as the flow of both fluids is governed by the Navier–Stokes equations, and the flow is coupled to dissolution and precipitation. Second, it is shown that, when letting the Cahn number approach zero, the limit of the two-scale phase-field model is an upscaled sharp-interface model. Finally, it is shown that the two processes, the upscaling and the sharp-interface limit, do commute. In other words, when starting with a diffuse-interface model at the pore scale, the order in which the Cahn number and the aspect ratio of the pore approach zero makes no difference, the result being in either case the upscaled counterpart of the sharp-interface model.

This paper is organized as follows. First, in § 2 a sharp-interface model for two fluid phases and one solid phase (including precipitation and dissolution) is presented. This model is approximated by the phase-field model proposed by Rohde & von Wolff (2021), which is discussed briefly in § 3. After bringing the phase-field model to a non-dimensional form in § 4, in § 5 we derive its upscaled counterpart by considering a thin strip geometry. The upscaled model still uses phase-field variables to locate the diffuse interfaces. In § 6 we identify the sharp-interface limit, that is the limit when letting the diffuse-interface width go to zero. Notably the upscaling and the sharp-interface limit commute. The numerical examples discussed in § 7 conclude the work.

2. The sharp-interface model

We start by presenting the sharp-interface model, which is then approximated by a phase-field model. For both the sharp-interface model and the phase-field model see Rohde & von Wolff (2021) for more details. We let $T > 0$ stand for the maximal time. For each $t \in [0, T]$, an $N$-dimensional domain $\varOmega$ ($N = 2$ or 3) is partitioned into three disjoint subdomains, $\varOmega _1(t), \varOmega _2(t)$ and $\varOmega _3(t)$. These are occupied by the two fluid phases and by the solid phase, respectively. The interface between the domain $\varOmega _i$ and $\varOmega _j$ is denoted by $\varGamma _{ij}$ ($i, j \in \{{1, 2, 3}\}, i \ne j$). Observe that these interfaces also depend on time.

With $t \in (0, T]$, in the fluid occupied subdomains $\varOmega _i(t), i\in \{{1,2}\}$ the model is governed by the incompressible Navier–Stokes equations

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {{\boldsymbol{v}}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\partial_t (\rho_i {{\boldsymbol{v}}}) + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \rho_i {{\boldsymbol{v}}} \otimes {{\boldsymbol{v}}} ) + \boldsymbol{\nabla} p = \boldsymbol{\nabla} \boldsymbol{\cdot} (2 \gamma_i \nabla^{s} {{\boldsymbol{v}}}), \end{gather}$$

where $\rho _i, \gamma _i$ denote the mass density, respectively viscosity of the fluid phase $i$, all assumed constant here; ${{\boldsymbol {v}}}$ and $p$ denote the fluid velocity and pressure in $\varOmega _i$, the index $i$ being skipped. The symmetrized strain (Jacobian) is given by $\nabla ^{s} {{\boldsymbol {v}}} = \frac {1}{2}(\boldsymbol {\nabla } {{\boldsymbol {v}}} + (\boldsymbol {\nabla } {{\boldsymbol {v}}})^{t})$.

At the interface $\varGamma _{12}(t)$ (separating $\varOmega _1(t)$ and $\varOmega _2(t)$) we assume that the velocity ${{\boldsymbol {v}}}$ is continuous and that the jump in the normal stress is only in the normal direction, and proportional to the curvature of the interface

(2.3)$$\begin{gather} {[}\![ {{{\boldsymbol{v}}}} ]\!] = 0, \end{gather}$$

(2.4)$$\begin{gather}{[}\![ {(p I - 2 \gamma \nabla^{s} {{\boldsymbol{v}}}) \boldsymbol{\cdot} {{\boldsymbol{n}}}} ]\!] = \sigma_{12} \kappa {{\boldsymbol{n}}}, \end{gather}$$

(2.5)$$\begin{gather}\nu = {{\boldsymbol{v}}} \boldsymbol{\cdot} {{\boldsymbol{n}}}. \end{gather}$$

Here, $[\![ {\cdot } ]\!]$ denotes the jump of a quantity over the interface, ${{\boldsymbol {n}}}$ the unit normal vector pointing outwards $\varOmega _1, \kappa$ the curvature of the interface and $\sigma _{12}$ the constant surface tension coefficient. Through the last condition, the normal velocity $\nu$ of the interface and the normal velocity of the fluids are equal.

The subdomain $\varOmega _3(t)$ is occupied by a mineral. We assume that the mineral phase is non-deforming and always connected to an outer boundary. Therefore, the mineral phase is stationary, and we do not need to solve for a velocity field ${{\boldsymbol {v}}}$ here. Note that this does not allow for small mineral grains that are transported by the fluid flow. The mineral phase is formed by the precipitation of two solute species present in fluid 1. The reverse process, in which the mineral can be dissolved and release solute in fluid 1 is also possible. In a simplified setting, assuming a constant electrical charge, it suffices to consider only one solute concentration in the model (see van Duijn & Knabner 1997), which is denoted by $c$. Here, we assume that solute is only present in fluid 1. Therefore, the solute transport is governed by the transport–diffusion equation in $\varOmega _1(t)$

(2.6)\begin{equation} \partial_t c + \boldsymbol{\nabla} \boldsymbol{\cdot} ({{\boldsymbol{v}}}c) - D {\rm \Delta} c = 0, \end{equation}

where $D$ is the constant diffusion coefficient.

The interface $\varGamma _{13}(t)$ is evolving due to precipitation and dissolution. At $\varGamma _{13}$ one has

(2.7)$$\begin{gather} \nu ={-} r(c) + \alpha \sigma_{13} \kappa, \end{gather}$$

(2.8)$$\begin{gather}D \boldsymbol{\nabla} c \boldsymbol{\cdot} {{\boldsymbol{n}}} = \nu (c^{{\ast}} - c). \end{gather}$$

The reaction rate $r(c)$ appearing in the former is generic and only depending on the solute concentration $c$. It accounts for dissolution and precipitation effects and is assumed increasing in $c$. In a more general situation we would expect the reaction rate to also depend on temperature. For the sake of simplicity we focus on the isothermal case here.

Remark 2.1 A simple reaction rate $r(c)$ can be constructed by assuming a constant dissolution rate $k_1$ and a quadratic mass action law, with rate $k_2$, for precipitation. With this the reaction rate is given by

(2.9)\begin{equation} r(c) ={-} k_1 + k_2 c^{2}. \end{equation}

The last term in (2.7), involving a constant parameter $\alpha \geq 0$ and the constant surface energy $\sigma _{13}$, allows for curvature effects in the evolution of $\varGamma _{13}$, e.g. an accelerated dissolution of solid tips peaking into the fluid phase. For $\alpha = 0$ no curvature effects enter the model.

Equation (2.8) is the Rankine–Hugoniot condition, ensuring the conservation of the solute species $c$. Here, $c^{\ast }$ is the (constant) molar concentration of the solute species, now as part of the mineral phase $\varOmega _3$. We refer to van Noorden (2009a) for the mathematical modelling details. For a specific application, we refer to von Wolff et al. (2021), where calcite precipitation is studied. There, $c$ denotes the molar concentration of inorganic carbon in the fluid, and $c^{\ast }$ is the molar density of calcium carbonate.

Equations (2.7) and (2.8) only hold at $\varGamma _{13}(t)$ and not at outer boundaries of $\varOmega$. That is, we do not allow for precipitation and dissolution at the outer boundaries of $\varOmega$.

At the fluid–fluid interface $\varGamma _{12}(t)$, a similar Rankine–Hugoniot condition is imposed

(2.10)\begin{equation} \boldsymbol{\nabla} c \boldsymbol{\cdot} {{\boldsymbol{n}}} = 0. \end{equation}

As before, ${{\boldsymbol {n}}}$ is the unit normal vector pointing outwards $\varOmega _1$. This condition ensures conservation of the solute $c$, since the concentration $c$ in fluid 2 is zero, and the normal velocity of the two fluids and of the interface are equal.

In contrast to $\varGamma _{13}$, no precipitation or dissolution are possible at the interface $\varGamma _{23}$ between $\varOmega _2$ and $\varOmega _3$. This is because we assume that fluid 2 does not contain any solute species. Therefore, the interface does not evolve, and its normal velocity is $\nu = 0$.

At the interfaces between a fluid and the mineral we impose the no-penetration condition

(2.11)\begin{equation} {{\boldsymbol{v}}} \boldsymbol{\cdot} {{\boldsymbol{n}}} = 0, \end{equation}

at $\varGamma _{12}$ and $\varGamma _{23}$.

Remark 2.2 Equation (2.11) ensures conservation of mass, under the following assumptions. First, we assume that the concentration of ions in the first fluid is small, so the mass density $\rho _1$ of this fluid phase is constant. Moreover, we assume that $\rho _1$ also equals the mass density $\rho _3$ of the solid phase. In the more general case of $\rho _1 \neq \rho _3$, a volume change associated with the reaction may appear, and one obtains the more general interface condition

(2.12)\begin{equation} -\rho_1 {{\boldsymbol{v}}} \boldsymbol{\cdot} {{\boldsymbol{n}}} = \nu [\![ {\rho} ]\!], \end{equation}

at $\varGamma _{12}$ and $\varGamma _{23}$. We refer to van Noorden (2009a) for more details on this aspect. For the sake of simplicity and because we assume that the normal velocity $\nu$ of the fluid–solid interfaces is small, we will present the simplified case of (2.11).

Finally, at the interfaces between a fluid and the mineral a Navier-slip condition (see Navier 1823) is assumed,

(2.13)\begin{equation} {{\boldsymbol{v}}} \boldsymbol{\cdot} {{{\tau}}} ={-}2 L_{slip} {{{\tau}}} (\nabla^{s} {{\boldsymbol{v}}}) {{\boldsymbol{n}}} , \end{equation}

at $\varGamma _{i3}$ ($i\in \{{1,2}\}$). Here, ${{{\tau }}} \in {{\mathbb {R}}}^{N}$ is any tangent vector to $\varGamma _{i3}$ (thus ${{{\tau }}} \perp {{\boldsymbol {n}}}$). The slip length $L_{slip} \geq 0$ is given by

(2.14)\begin{equation} L_{slip} = \gamma_1 \sqrt{\frac{2}{ \rho_3 d_0 \gamma_3}}. \end{equation}

Here, $\rho _3$ is the density of the solid phase, and the constant $d_0$ will be determined by choices in the phase-field model. As explained in Rohde & von Wolff (2021), $\gamma _3$ is not the viscosity of the solid phase, but can be chosen instead to archive a given slip length $L_{slip}$ in the relation (2.14). While the model allows for a positive slip length to account for additional properties of the fluid–solid interface, e.g. surface roughness, one can also choose $\gamma _3$ large enough to get $L_{slip} \approx 0$.

As a last step we consider points where the fluid–fluid interface $\varGamma _{12}$ intersects the solid boundary $\varGamma _{13} \cup \varGamma _{23}$. At these contact points all bulk domains $\varOmega _1, \varOmega _2, \varOmega _3$ meet. We assume that the set of contact points consist of distinct points when in two dimensions and of distinct lines in the three-dimensional case. We focus first on the two-dimensional case.

Given the constant surface energies $\sigma _{12}, \sigma _{13}, \sigma _{23} > 0$ we impose for the contact angle the condition

(2.15)\begin{equation} \frac{\sin(\beta_1)}{\sigma_{23}} = \frac{\sin(\beta_2)}{\sigma_{13}} = \frac{\sin(\beta_3)}{\sigma_{12}}, \end{equation}

with $\beta _i$ being the contact angle of phase $\varOmega _i$ at the contact point. Together with the condition $\beta _1 + \beta _2 + \beta _3 = 2{\rm \pi}$ this uniquely determines the contact angles $\beta _i$. In the three-dimensional case, the same condition (2.15) is imposed on the plane perpendicular to the contact line.

3. The phase-field model

The sharp-interface model in § 2 involves free boundaries, which makes it difficult from both the analysis and numerical points of view. Relying on the idea to approximate the characteristic functions of each of the phases by smooth phase indicators (Caginalp & Fife 1988), phase-field models are convenient alternatives. For the specific problem considered here, a phase-field model called $\delta$-$2f1s$-model was introduced in Rohde & von Wolff (2021); here, we present it briefly for completeness. We refer to Rohde & von Wolff (2021) for more details on the derivation and the properties of the model, including the derivation of the sharp-interface limit.

3.1. Preliminaries

The $\delta$-$2f1s$-model introduces three phase-field variables $\phi _1, \phi _2, \phi _3$ that represent the volume fraction of the two fluid phases and of the solid phase, respectively. Thus, $\phi _i$ approximates the indicator function of $\varOmega _i$ appearing in the sharp-interface model in § 2. The phase-field variables ${{\boldsymbol {\varPhi }}} = (\phi _1, \phi _2, \phi _3)^{t}$ are smooth and defined on the entire domain $\varOmega$. In the sharp-interface model, the transition from one phase to another is across an interface. In the phase-field model, this interface is replaced by a diffuse transition zone from one phase to another, where the gradients of the corresponding phase-field variables are high. A ternary Cahn–Hilliard equation governs the evolution of ${{\boldsymbol {\varPhi }}}$, and is coupled with a Navier–Stokes equation for fluid flow, and a reaction–transport–diffusion equation for dissolved ion concentration $c$.

The $\delta$-$2f1s$-model additionally introduces a small regularization parameter $\delta > 0$. Since no maximum principle holds for the Cahn–Hilliard equation, $\delta$ is used to ensure the positivity of the volume fractions. Also, the double-well potential

(3.1)\begin{align} W_{{dw}}(\phi) &= 450 \phi^{4} (1-\phi)^{4} + \delta \ell\left(\frac{\phi}{\delta} \right) + \delta \ell \left( \frac{1-\phi}{\delta} \right),\nonumber\\ &\quad \text{with } \ell(x) = \begin{cases} \dfrac{x^{2}}{1+x} & x \in ({-}1,0), \\ 0 & x \geq 0, \end{cases} \end{align}

is employed. Observe that $W_{dw}$ has two minima at $0$ and $1$, and becomes unbounded at $-\delta$ and $1+\delta$. With this, we define the triple-well potential

(3.2)\begin{equation} W({{\boldsymbol{\varPhi}}}) := W_0(P {{\boldsymbol{\varPhi}}}), \quad \text{where } W_0({{\boldsymbol{\varPhi}}}) = \sum_{i=1}^{3} \varSigma_i W_{dw}(\phi_i). \end{equation}

Here, $\varSigma _i > 0$ are surface energy coefficients, and $P$ is the projection of ${{\mathbb {R}}}^{3}$ onto the plane $\sum _i \phi _i = 1$, given by

(3.3a,b)\begin{equation} P{{\boldsymbol{\varPhi}}} = {{\boldsymbol{\varPhi}}} + \varSigma_T(1-\phi_1-\phi_2-\phi_3) \begin{pmatrix} \varSigma_1^{{-}1}\\ \varSigma_2^{{-}1}\\ \varSigma_3^{{-}1} \end{pmatrix} ,\quad \frac{1}{\varSigma_T} = \frac{1}{\varSigma_1}+\frac{1}{\varSigma_2}+ \frac{1}{\varSigma_3}. \end{equation}

As shown in Rohde & von Wolff (2021), this construction ensures that the volume fractions sum to one, i.e. $\sum _{i=1}^{3} \phi _i = 1$, provided the initial data have this property. Furthermore, Rohde & von Wolff (2021) uses an energy argument and the unboundedness of the potential to show that $-\delta < \phi _i < 1 + \delta$ ($i = 1, 2, 3$).

Next, we define the total fluid volume fraction $\tilde \phi _f$ and ion-dissolving fluid fraction $\phi _c$ as

(3.4)$$\begin{gather} \tilde \phi_f := \phi_1+\phi_2 + 2\delta \phi_3, \end{gather}$$

(3.5)$$\begin{gather}\phi_c := \phi_1, \end{gather}$$

(3.6)$$\begin{gather}\tilde \phi_c := \phi_1 + \delta. \end{gather}$$

Here, the tilde denotes a modification using the small parameter $\delta$, to ensure that the respective variables are positive. Using the (constant) fluid densities $\rho _i$ and viscosities $\gamma _i$ the total fluid density $\rho _f$ and viscosity $\tilde \gamma$ become

(3.7)$$\begin{gather} \rho_f({{\boldsymbol{\varPhi}}}) := \rho_1 \phi_1 + \rho_2 \phi_2, \end{gather}$$

(3.8)$$\begin{gather}\tilde \rho_f({{\boldsymbol{\varPhi}}}) := \rho_1 \phi_1 + \rho_2 \phi_2 + (\rho_1+\rho_2) \delta, \end{gather}$$

(3.9)$$\begin{gather}\tilde \gamma({{\boldsymbol{\varPhi}}}) := (\phi_1 \gamma_1^{{-}1}+ \phi_2 \gamma_2^{{-}1} + \phi_3 \gamma_3^{{-}1} + (\gamma_1^{{-}1}+\gamma_2^{{-}1}+\gamma_3^{{-}1})\delta)^{{-}1}. \end{gather}$$

3.2. The $\delta$-$2f1s$-model

We now present the $\delta$-$2f1s$-model. All equations are defined in $(0, T] \times \varOmega$. The flow is governed by the Navier–Stokes equations and involves the fluid fraction $\tilde \phi _f$,

(3.10)\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} (\tilde \phi_f {{\boldsymbol{v}}}) = 0, \end{align}

(3.11)\begin{align} \partial_t (\tilde\rho_f {{\boldsymbol{v}}}) + \boldsymbol{\nabla} \boldsymbol{\cdot} ((\rho_f {{\boldsymbol{v}}} + \rho_1 {{\boldsymbol{J}}}_1 + \rho_2 {{\boldsymbol{J}}}_2) \otimes {{\boldsymbol{v}}}) &={-} \tilde\phi_f \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} (2 \tilde \gamma({{\boldsymbol{\varPhi}}}) \nabla^{s} {{\boldsymbol{v}}}) \nonumber\\ &\quad - \rho_3 d(\tilde\phi_f) {{\boldsymbol{v}}} + \tilde {{\boldsymbol{S}}} + \frac 12 \rho_1 {{\boldsymbol{v}}}R. \end{align}

This is coupled with the transport–diffusion–reaction equation for the ion concentration

(3.12)\begin{equation} \partial_t (\tilde\phi_c c) + \boldsymbol{\nabla} \boldsymbol{\cdot} ((\phi_c {{\boldsymbol{v}}} + {{\boldsymbol{J}}}_1) c) = D \boldsymbol{\nabla} \boldsymbol{\cdot} ( \tilde\phi_c \boldsymbol{\nabla} c) + c^{{\ast}} R. \end{equation}

The phase-field variables $\phi _1, \phi _2, \phi _3$ satisfy the Cahn–Hilliard equations

(3.13)$$\begin{gather} \partial_t \phi_1 + \boldsymbol{\nabla} \boldsymbol{\cdot} (\phi_1 {{\boldsymbol{v}}} + {{\boldsymbol{J}}}_1) = R, \end{gather}$$

(3.14)$$\begin{gather}\partial_t \phi_2 + \boldsymbol{\nabla} \boldsymbol{\cdot} (\phi_2 {{\boldsymbol{v}}} + {{\boldsymbol{J}}}_2) = 0, \end{gather}$$

(3.15)$$\begin{gather}\partial_t \phi_3 + \boldsymbol{\nabla} \boldsymbol{\cdot} (2\delta \phi_3 {{\boldsymbol{v}}} + {{\boldsymbol{J}}}_3) ={-}R, \end{gather}$$

(3.16)$$\begin{gather}{{\boldsymbol{J}}}_i ={-} \frac{{{\varepsilon}}M}{\varSigma_i} \boldsymbol{\nabla} \mu_i, \quad i\in\{{1,2,3}\}, \end{gather}$$

(3.17)$$\begin{gather}\mu_i = \frac{\partial_{\phi_i}W({{\boldsymbol{\varPhi}}})}{{{\varepsilon}}} - {{\varepsilon}}\varSigma_i {\rm \Delta} \phi_i, \quad i\in\{{1,2,3}\}. \end{gather}$$

Compared with the common Navier–Stokes equations, some modifications appear in (3.11). The fluid density $\tilde \rho _f({{\boldsymbol {\varPhi }}})$ introduces a strong coupling between the Navier–Stokes equations and the Cahn–Hilliard equations. All terms except the advection term use the modified quantities $\tilde \phi _f, \tilde \rho _f$ and $\tilde \gamma$. Additional flux terms $\rho _i {{\boldsymbol {J}}}_i \otimes {{\boldsymbol {v}}}$ are introduced to account for momentum fluxes due to the Cahn–Hilliard evolution. Secondly, the dissipative term $-\rho _3 d(\tilde \phi _f) {{\boldsymbol {v}}}$ is added. Here, $d$ is a decreasing function such that $d(0) = d_0 > 0$ and $d(1)=0$, for example $d(\tilde \phi _f) = d_0 (1-\tilde \phi _f)^{2}$. The term $d(\tilde \phi _f)$ is therefore active in the solid phase and guarantees that ${{\boldsymbol {v}}}$ remains small there. It also influences the slip length $L_{slip}$. Lastly, the surface tension term $\tilde {{\boldsymbol {S}}}$ is given by

(3.18)\begin{equation} \tilde {{\boldsymbol{S}}} ={-}\mu_2 \tilde \phi_f\boldsymbol{\nabla}\left(\frac{\phi_1}{\tilde\phi_f}\right) -\mu_1 \tilde\phi_f\boldsymbol{\nabla}\left(\frac{\phi_2}{\tilde\phi_f}\right) - 2\delta \phi_3 \boldsymbol{\nabla} (\mu_3-\mu_1-\mu_2) . \end{equation}

The reaction term $R$ modelling precipitation and dissolution of ions is given by

(3.19)\begin{equation} R ={-}q({{\boldsymbol{\varPhi}}}) (r(c)+ \tilde \alpha \mu_1-\tilde \alpha \mu_3). \end{equation}

Here, $r(c)$ is the increasing reaction rate used in the sharp-interface description (2.7). Additionally, the precipitation process can depend on curvature effects through surface effects that are similar to surface diffusion, and are encountered if $\alpha > 0$. Again, the tilde denotes a modification of $\alpha$, that is $\tilde \alpha = \alpha + \delta$. Finally, to concentrate the reaction inside the diffuse interface region between fluid phase 1 and the solid phase, which is equivalent to the assumption made in the sharp-interface model, the non-dimensional term $q({{\boldsymbol {\varPhi }}}) = 30 \phi _1^{2} \phi _3^{2}$ is used. Observe that $q$ dominates wherever neither $\phi _1$ nor $\phi _2$ are close to 0, which is precisely the envisaged location for the fluid 1–mineral interface.

4. Non-dimensionalization

We proceed by bringing the $\delta$-$2f1s$-model (3.10)–(3.17) to a non-dimensional form, and derived an upscaled counterpart of it by employing asymptotic expansion and averaging techniques. We consider a simplified geometric setting. We start by introducing a thin strip having length $L$ and width $\ell \ll L$, as shown in figure 1.

Figure 1. Setting of the thin strip: the strip with length $L$ and width $\ell$ consists of solid walls (red, ${{\boldsymbol {\varPhi }}} \approx (0,0,1)^{t}$) and fluid phases (light blue, dark blue). The diffuse-interface width ${{\varepsilon }}$ is smaller than $\ell$.

With a chosen domain width $\ell _\varOmega > \ell$, the domain $\varOmega = [0,L] \times [-\ell _\varOmega /2, \ell _\varOmega /2]$ includes the thin strip mentioned above, which is identified as $[0, L] \times [-\ell / 2, \ell / 2]$. The region outside the strip is occupied by the mineral, so ${{\boldsymbol {\varPhi }}} \approx (0,0,1)^{t}$ there. The diffuse interfaces are located in regions of width ${{\varepsilon }}$. We assume here that the diffuse-interface regions remain clearly separated inside the thin strip, hence ${{\varepsilon }} \ll \ell$.

Three length scales can be identified, $L\gg \ell \gg {{\varepsilon }}$. These are related through the aspect ratio $\beta = \ell /L$ and the Cahn number $Cn = {{\varepsilon }}/L$, both assumed small. Observe that, in fact, $Cn \ll \beta \ll 1$.

The reference quantities used in the non-dimensionalization procedure are listed in table 1. Non-dimensional values are then identified by a hat. Note that we relate only few reference values directly to each other. In particular, we do relate reference values when we want to emphasize an explicit dependence on ${y}_{{ref}}$, as seen for ${p}_{{ref}}, {d}_{{ref}}$ and ${\mu }_{{ref}}$. The choices are motivated as follows. To obtain an upscaled macroscopic velocity of order ${v}_{{ref}} = {x}_{{ref}} / {t}_{{ref}}$, the pressure drop in the thin strip has to scale with $1/({y}_{{ref}})^{2}$. Also, the slip length $L_{slip}$ is supposed to be of order $\ell$ and not $L$, which is achieved by a momentum dissipation scaling $1/({y}_{{ref}})^{2}$.

Table 1. Variables, reference values and non-dimensional quantities for the non-dimensionalization.

We rewrite the Cahn number introduced above in terms of reference quantities, and define other dimensionless numbers that are used below to relate the reference quantities: the Reynolds number, Capillary number, Damköhler number and Péclet numbers of the Cahn–Hilliard (CH) evolution and ion concentration

(4.1)\begin{equation} \left.\begin{array}{c@{}} Re = \dfrac{ {\rho}_{{ref}} {v}_{{ref}} {x}_{{ref}}}{{\gamma}_{{ref}}}, \quad Ca = \dfrac{ {\gamma}_{{ref}} {v}_{{ref}}}{{\varSigma}_{{ref}}}, \quad Cn = \dfrac{ {{\varepsilon_{ref}}}}{{x}_{{ref}}},\\ Da = \dfrac{ {r}_{{ref}} {x}_{{ref}}}{{v}_{{ref}}}, \quad Pe_{CH} = \dfrac{ {v}_{{ref}} {x}_{{ref}}}{ M_{ref}}, \quad Pe_{c} = \dfrac{ {v}_{{ref}} {x}_{{ref}}}{D_{ref}}. \end{array}\right\} \end{equation}

Clearly, the non-dimensionalization also affects the spatial and temporal derivatives, namely

(4.2a,b)\begin{equation} \hat {\boldsymbol{\nabla}} = {x}_{{ref}} \boldsymbol{\nabla} \quad \text{and} \quad \partial_{\hat t} = {t}_{{ref}} \partial_t. \end{equation}

We now can insert the non-dimensional variables of table 1, the non-dimensional numbers (4.1) and the non-dimensional operators in (4.2a,b) into the $\delta$-$2f1s$-model (3.10)–(3.17). The non-dimensional equations become

(4.3)\begin{gather} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\tilde \phi_f \hat{{\boldsymbol{v}}}) = 0,\end{gather}

(4.4)\begin{gather} \partial_{\hat t} (\hat{\tilde\rho}_f \hat{{\boldsymbol{v}}}) + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\hat{\rho}_f \hat{{\boldsymbol{v}}} \otimes \hat{{\boldsymbol{v}}}) + \frac{Cn}{\beta Pe_{CH}} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} ((\hat{\rho}_1 \hat{{\boldsymbol{J}}}_1 + \hat{\rho}_2 \hat{{\boldsymbol{J}}}_2) \otimes \hat{{\boldsymbol{v}}}) \nonumber\\ \hspace{-5.5pc}={-} \frac{1}{\beta^{2} Re}\tilde\phi_f \hat{\boldsymbol{\nabla}} \hat{p} + \frac{1}{Re} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (2 \hat{\tilde \gamma}({{\boldsymbol{\varPhi}}}) \hat{\nabla}^{s} \hat{{\boldsymbol{v}}}) \nonumber\\ - \frac{1}{\beta^{2} Re}\hat \rho_3 \hat d(\tilde\phi_f) \hat{{\boldsymbol{v}}} + \frac{1}{\beta Re} \frac{1}{Ca} \hat{\tilde {{\boldsymbol{S}}}} + Da \frac 12 \hat{\rho}_1 \hat{{\boldsymbol{v}}} \hat{R}, \end{gather}

for the flow,

(4.5)\begin{equation} \partial_{\hat t} (\tilde\phi_c \hat c) + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\phi_c \hat{{\boldsymbol{v}}} \hat c) + \frac{Cn}{\beta Pe_{CH}} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot}(\hat{{\boldsymbol{J}}}_1 \hat c) = \frac{1}{Pe_c} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} ( \tilde\phi_c \hat{\boldsymbol{\nabla}} \hat c) + Da \hat R,\end{equation}

for the ion transport–diffusion–reaction, while for the Cahn–Hilliard evolution one gets

(4.6)$$\begin{gather} \partial_{\hat t} \phi_1 + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\phi_1 \hat{{\boldsymbol{v}}}) + \frac{Cn}{\beta Pe_{CH}} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \hat{{\boldsymbol{J}}}_1 = Da \hat R, \end{gather}$$

(4.7)$$\begin{gather}\partial_{\hat t} \phi_2 + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\phi_2 \hat{{\boldsymbol{v}}}) + \frac{Cn}{\beta Pe_{CH}} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \hat{{\boldsymbol{J}}}_2 = 0, \end{gather}$$

(4.8)$$\begin{gather}\partial_{\hat t} \phi_3 + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} (2\delta \phi_3 \hat{{\boldsymbol{v}}}) + \frac{Cn}{\beta Pe_{CH}} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \hat{{\boldsymbol{J}}}_3 ={-}Da \hat R, tag{4.8}\end{gather}$$

(4.9)$$\begin{gather}\hat{{\boldsymbol{J}}}_i ={-} \frac{1}{\hat \varSigma_i} \hat{\boldsymbol{\nabla}} \hat \mu_i, \quad i\in\{{1,2,3}\}, \end{gather}$$

(4.10)$$\begin{gather}\frac{\hat \mu_i}{\beta } = \frac{\partial_{\phi_i}\hat W({{\boldsymbol{\varPhi}}})}{Cn} - Cn \hat \varSigma_i \hat \nabla^{2} \phi_i, \quad i\in\{{1,2,3}\}. \end{gather}$$

All equations are defined in the dimensionless time–space domain $(0, 1] \times \hat {\varOmega }$, where $\hat \varOmega = [0,1] \times [- \hat \ell _\varOmega /2, \hat \ell _\varOmega /2]$. The surface tension and reaction are given as

(4.11)$$\begin{gather} \hat{\tilde {{\boldsymbol{S}}}}={-}\hat \mu_2 \tilde \phi_f\hat{\boldsymbol{\nabla}}\left(\frac{\phi_1}{\tilde\phi_f}\right) -\hat \mu_1 \tilde\phi_f\hat{\boldsymbol{\nabla}}\left(\frac{\phi_2}{\tilde\phi_f}\right) - 2\delta \phi_3 \hat{\boldsymbol{\nabla}} (\hat \mu_3-\hat \mu_1-\hat \mu_2), \end{gather}$$

(4.12)$$\begin{gather}\hat R ={-}q({{\boldsymbol{\varPhi}}}) (\hat r(\hat c) + \hat{\tilde \alpha} \hat \mu_1 - \hat{\tilde \alpha} \hat \mu_3). \end{gather}$$

From here on, we will only work with the non-dimensional model and therefore the hats are left out in the notation.

5. Upscaling in a thin strip

We now proceed by deriving the upscaled model, obtained when passing to the limit $\beta \to 0$. This means that the thin strip reduces to a one-dimensional object, as its width is vanishing compared with its length.

We introduce new coordinates $(x,y)$ such that ${{\boldsymbol {x}}} = (x, \beta y)$. In the thin strip we expect all variables to vary in the longitudinal direction ${{\boldsymbol {e}}}_x$ on the length scale $L = {x}_{{ref}}$ and in the transverse direction ${{\boldsymbol {e}}}_y$ on the length scale $\ell = {y}_{{ref}} = \beta {x}_{{ref}}$. In particular, this will result in $\boldsymbol {\nabla } = {{\boldsymbol {e}}}_x \partial _x + \beta ^{-1} {{\boldsymbol {e}}}_y \partial _y$.

The non-dimensional domain is given by $\varOmega = [0,1]\times [-\ell _\varOmega /2, \ell _\varOmega /2]$ (recall that we dropped the hats in the notation) and we choose for the upscaling the boundary conditions at $y = \pm \ell _\varOmega /2$ as

(5.1)$$\begin{gather} \partial_y {{\boldsymbol{\varPhi}}}(t,x,\pm \ell_\varOmega/2) = 0, \end{gather}$$

(5.2)$$\begin{gather}\partial_y \mu(t,x,\pm \ell_\varOmega/2) = 0, \end{gather}$$

(5.3)$$\begin{gather}\partial_y c(t,x,\pm \ell_\varOmega/2) = 0, \end{gather}$$

(5.4)$$\begin{gather}{{\boldsymbol{v}}}(t,x,\pm \ell_\varOmega/2) = 0. \end{gather}$$

5.1. Scaling of non-dimensional numbers

The upscaled model will also depend on the scaling of the dimensionless numbers (4.1) with respect to $\beta$. We consider the following behaviour of these numbers with respect to $\beta$

(5.5)$$\begin{gather} Re = \overline{Re}, \end{gather}$$

(5.6)$$\begin{gather}Ca = \overline{Ca}, \end{gather}$$

(5.7)$$\begin{gather}Cn = \beta \bar {{\varepsilon}}, \end{gather}$$

(5.8)$$\begin{gather}Da = \overline{Da} / \bar {{\varepsilon}}, \end{gather}$$

(5.9)$$\begin{gather}Pe_{CH} = 1 / (\beta^{2} \bar M), \end{gather}$$

(5.10)$$\begin{gather}Pe_{c} = \overline{Pe_c}, \end{gather}$$

where $\overline {Re}, \overline {Ca}, \bar {{\varepsilon }}, \overline {Da}, \bar M, \overline {Pe_c}$ are constants independent of $\beta$. In detail, these choices are motivated as follows.

• • The moderate Reynolds number (5.5) leads to a parabolic flow profile in the thin strip, we expect laminar flow.

• • As the curvature of the fluid–fluid interface is of order $O(\beta )$, choosing a moderate capillary number $Ca$ in (5.6) leads to the same pressure in both fluids, thus the capillary pressure becomes 0 (for sharp-interface models see also Sharmin et al. 2020; Lunowa et al. 2021). Note that this is a major difference to the three-dimensional case (see e.g. Mikelić 2009) where we expect a curvature of $O(\beta ^{-1})$ leading to a non-zero capillary pressure.

• • The scaling of the Cahn number $Cn$ in (5.7) can be reformulated to $\bar {{\varepsilon }} = {{\varepsilon }} / \ell$. Therefore, the interface width ${{\varepsilon }}$ scales with the width of the thin strip, $\ell$. At the same time, the diffuse-interface regions are assumed to be localized inside the thin strip, therefore we require ${{\varepsilon }} \ll \ell$. This translates into a fixed, small $\bar {{\varepsilon }}$, i.e. $\bar {{\varepsilon }} \ll 1$. In the numerical experiments presented in § 7 we choose $\bar {{\varepsilon }} = 0.03$.

• • We consider a moderate Damköhler number (5.8). In the sharp-interface model, this would ensure that the interfaces move with moderate velocity inside the thin strip, proportional to $\ell /T$. In the diffuse-interface model, the reaction is only active in the diffuse-interface region, which has an area scaling with ${{\varepsilon }}$. Therefore, $Da$ is divided by $\bar {{\varepsilon }}$, and expect to have fluid–solid or fluid–fluid interfaces evolving over the length scale $\ell$. A dominating Damk’ ohler regime like $Da = O(\beta ^{-1})$ would instead lead to equilibrium-type reactions in the upscaled model, but the evolution of the interfaces should remain moderate. This can be achieved by assuming that the molar density of the species in the precipitate is sufficiently high to compensate the fast reaction kinetics.

• • The high Péclet number (5.9) for the phase field assures that the evolution of the phase field remains within the transverse length scale $\ell$ in an $O(1)$ time scale.

• • The moderate Péclet number of the ion diffusion (5.10) will result in a macroscopic diffusion of ions, while the ion distribution in the transverse direction equilibrates faster than the $O(1)$ time scale.

Lastly, the small, non-dimensional number $\delta >0$ appears in the $\delta$-$2f1s$-model. It is used as a regularization parameter, to ensure the positivity of volume fractions, density and viscosity. Here, we assume that $\delta$ is constant and independent of $\beta$.

5.2. Asymptotic expansions

We assume that we can write solutions to the non-dimensional $\delta$-$2f1s$-model (4.3)–(4.10) in terms of an asymptotic expansion in $\beta$ of ${{\boldsymbol {\varPhi }}}, {{\boldsymbol {v}}}, p, c, \mu _1, \mu _2, \mu _3$. To be precise, we assume expansions of the form

(5.11)\begin{equation} {{\boldsymbol{\varPhi}}}(t,{{\boldsymbol{x}}}) = {{\boldsymbol{\varPhi}}}_0(t,x,y) + \beta {{\boldsymbol{\varPhi}}}_1(t,x,y) + \beta^{2} {{\boldsymbol{\varPhi}}}_2(t,x,y) + \cdots, \end{equation}

where ${{\boldsymbol {\varPhi }}}_k, k \in {{\mathbb {N}}}_0$ does not depend on $\beta$.

Inserting these asymptotic expansions into the non-dimensional $\delta$-$2f1s$-model we group by powers of $\beta$. As the calculations are lengthy, we show them in Appendix A.

Remark 5.1 Note that the asymptotic expansions are written depending on the new coordinates $x$ and $y$. This means that in the ${{\boldsymbol {e}}}_x$ direction variables cannot vary on the (non-dimensional) length scale $\beta$, because a non-trivial function $f(x/\beta )$ cannot be expanded in the form $f(x/\beta ) = f_0(x) + \beta f_1(x) + \cdots$. In particular, this implies that there are no phase-field interfaces possible perpendicular to the thin strip, as they would change the value of ${{\boldsymbol {\varPhi }}}$ over the length $Cn = \beta \bar {{\varepsilon }}$. We will discuss in § 7.2 a numerical example that violates this assumption.

The assumption is also violated for triple points, where all three phases meet, and for points where interfaces meet the boundary of $\varOmega$ at $y = \pm \ell _\varOmega /2$. Therefore, $\ell _\varOmega$ has to be chosen big enough, such that the width of the thin strip does not reach $\ell _\varOmega$.

We will present in § 8 some ideas to handle cases where the assumption of slow variation in the ${{\boldsymbol {e}}}_x$ direction is violated.

5.3. The upscaled $\delta$-$2f1s$-model

Let us summarize the results of the upscaling done in detail in Appendix A. Except for ${{\boldsymbol {v}}}$ we will only need the leading-order term of each unknown, and will therefore drop the subscript $0$. We will call the model (5.12)–(5.27) the upscaled $\delta$-$2f1s$-model.

From (A6) and (A21) we have the macroscopic continuity equation for the total flux $Q_f$ and the Darcy equation for the pressure $p$, and the macroscopic transport–diffusion–reaction equation for the ion concentration $c$ (A18)

(5.12)$$\begin{gather} \partial_x Q_f = 0, \end{gather}$$

(5.13)$$\begin{gather}Q_f ={-}K_f \partial_x p, \end{gather}$$

(5.14)$$\begin{gather}\frac{{{\rm d}}}{{{\rm d}}t} ( \tilde \phi_{c,{total}} c ) + \partial_x ( ({-}K_c \partial_x p) c ) = \frac{1}{\overline{Pe_c}} \partial_x ( \tilde \phi_{c,{total}} \partial_x c) + \frac{\overline{Da}}{\bar {{\varepsilon}} } R_{total}. \end{gather}$$

These equations are macroscopic in the sense that the unknowns $Q_f, p$ and $c$ depend only on $x$ and $t$, but not on $y$. The parameters in these equations are upscaled quantities, depending on the exact distribution of the phases in the $y$ direction

(5.15)$$\begin{gather} \tilde \phi_{c,{total}} = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \tilde \phi_{c}\, {{\rm d}}y, \end{gather}$$

(5.16)$$\begin{gather}K_f(t,x) = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \tilde \phi_{f} w\, {{\rm d}}y, \end{gather}$$

(5.17)$$\begin{gather}K_c(t,x) = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \tilde \phi_{c} w\, {{\rm d}}y, \end{gather}$$

(5.18)$$\begin{gather}R_{total} = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} R\, {{\rm d}}y. \end{gather}$$

For the phase-field parameters we still have to solve the fully coupled two-dimensional problem (A7), (A9), (A10), (A11), that is

(5.19)$$\begin{gather} \partial_{t} \phi_{1} + \partial_x (\phi_{1} {{\boldsymbol{v}}}_0^ {(1)}) + \partial_y (\phi_{1} {{\boldsymbol{v}}}_1^ {(2)}) - \frac{\bar {{\varepsilon}} \bar M}{\varSigma_1} \partial_y^{2} \mu_{1} = \frac{\overline{Da}}{\bar {{\varepsilon}} }R, \end{gather}$$

(5.20)$$\begin{gather}\partial_{t} \phi_{2} + \partial_x (\phi_{2} {{\boldsymbol{v}}}_0^ {(1)}) + \partial_y (\phi_{2} {{\boldsymbol{v}}}_1^ {(2)}) - \frac{\bar {{\varepsilon}} \bar M}{\varSigma_2} \partial_y^{2} \mu_{2} = 0, \end{gather}$$

(5.21)$$\begin{gather}\partial_{t} \phi_{3} + \partial_x (2\delta \phi_{3} {{\boldsymbol{v}}}_0^ {(1)}) + \partial_y (2\delta \phi_{3} {{\boldsymbol{v}}}_1^ {(2)}) - \frac{\bar {{\varepsilon}} \bar M}{\varSigma_3} \partial_y^{2} \mu_{3} ={-}\frac{\overline{Da}}{\bar {{\varepsilon}}}R, \end{gather}$$

(5.22)$$\begin{gather}\mu_{i} = \frac{\partial_{\phi_i} W({{\boldsymbol{\varPhi}}})}{\bar {{\varepsilon}}} - \bar {{\varepsilon}} \varSigma_i \partial_y^{2} \phi_{i}, \quad i\in\{{1,2,3}\}, \end{gather}$$

with the reaction term

(5.23)\begin{equation} R ={-} q({{\boldsymbol{\varPhi}}})(r(c) + \tilde \alpha \mu_{1} - \tilde \alpha \mu_{3}). \end{equation}

Note that, in contrast to the non-dimensional model (4.3)–(4.10), the Cahn–Hilliard evolution acts only in the ${{\boldsymbol {e}}}_y$ direction. The only term acting in the ${{\boldsymbol {e}}}_x$ direction is the transport of the fluid phases. This will enable us in § 7.1 to develop a numerical model that uses explicit upwinding for the fluid transport and can therefore decouple cell problems for different values of $x$.

For the flow it suffices to solve the cell problem (A22), (A23)

(5.24)$$\begin{gather} \rho_3 d(\tilde\phi_{f}) w - \partial_y (\tilde \gamma({{\boldsymbol{\varPhi}}}) \partial_y w) = \tilde\phi_{f}, \end{gather}$$

(5.25)$$\begin{gather}\lim_{y\to\pm {\ell_\varOmega/2}} w = 0, \end{gather}$$

and recover the flow ${{\boldsymbol {v}}}^{(1)}_0, {{\boldsymbol {v}}}^{(2)}_1$ by (A21) and (A5)

(5.26)$$\begin{gather} {{\boldsymbol{v}}}^{(1)}_0 ={-}w \partial_x p, \end{gather}$$

(5.27)$$\begin{gather}\partial_x (\tilde \phi_{f} {{\boldsymbol{v}}}^{(1)}_0) + \partial_y (\tilde \phi_{f} {{\boldsymbol{v}}}^{(2)}_1) = 0. \end{gather}$$

Note that, while the equations for the flow (5.24)–(5.27) do not explicitly depend on time, they depend on the phase-field parameters ${{\boldsymbol {\varPhi }}}$, which can change in time.

6. Sharp-interface limit of the upscaled $\delta$-$2f1s$-model

In the previous section we have investigated the scale separation $\beta = \ell / L \to 0$. A different limit process that is commonly investigated for phase-field models is the sharp-interface limit ${{\varepsilon }} \to 0$. In Rohde & von Wolff (2021) this limit is analysed for the $\delta$-$2f1s$-model (3.10)–(3.17), resulting in the sharp-interface evolution described in § 2.

Because the upscaled $\delta$-$2f1s$-model (5.12)–(5.27) still contains a Cahn–Hilliard evolution, depending on the small number $\bar {{\varepsilon }} = {{\varepsilon }} / \ell$, we can investigate the sharp-interface limit $\bar {{\varepsilon }} \to 0$ of the upscaled $\delta$-$2f1s$-model. This means that we are interested in the limit process of vanishing diffuse-interface width ${{\varepsilon }}$ compared with the width $\ell$ of the thin strip. In the following, we will use matched asymptotic expansions to analyse this limit, the argumentation is mostly analogous to Rohde & von Wolff (2021).

6.1. Assumptions and scaling of non-dimensional numbers

To derive the sharp-interface limit $\bar {{\varepsilon }} \to 0$, we assume that $\overline {Pe_c}, \overline {Da}, \bar {M}$ are constant and independent of $\bar {{\varepsilon }}$. This choice of scaling allows for a reasonable limit process, with physical properties independent of the diffuse interface width.

The scaling $\delta = \bar {{\varepsilon }}$ is important. The regularization parameter $\delta$ is introduced in the $\delta$-$2f1s$-model to ensure the positivity of e.g. the density $\tilde \rho _f({{\boldsymbol {\varPhi }}})$ in (3.8). This $\delta$-regularization is not necessary for the sharp-interface formulation, and the choice $\delta = \bar {{\varepsilon }}$ leads to $\delta$ vanishing in the sharp-interface limit.

As a basic assumption we expect to have solutions that form bulk phases, characterized by nearly constant ${{\boldsymbol {\varPhi }}}$, and interfaces, characterized by a large gradient of ${{\boldsymbol {\varPhi }}}$. We also assume that $\mu _i, i\in \{{1,2,3}\}$ is of $O(1)$, not of $O(\bar {{\varepsilon }}^{-1})$, as (5.22) would suggest. For a discussion of why this assumption is reasonable on a $O(1)$ time scale, see Pego & Penrose (1989).

We also assume that in an interface between phase ${{\boldsymbol {\varPhi }}} = {{\boldsymbol {e}}}_i$ and ${{\boldsymbol {\varPhi }}} = {{\boldsymbol {e}}}_j$ the third phase is not present. This assumption is reasonable because with our constructions of $W$ (3.2) minimizers of the Ginzburg–Landau energy $W({{\boldsymbol {\varPhi }}}) + \sum _i \varSigma _i {\rm \Delta} \phi _i$ that connect ${{\boldsymbol {\varPhi }}} = {{\boldsymbol {e}}}_i$ and ${{\boldsymbol {\varPhi }}} = {{\boldsymbol {e}}}_j$ satisfy $\phi _k = 0, k\in \{{1,2,3}\}\setminus \{{i,j}\}$.

As the calculations for the sharp-interface limit are lengthy, we present them in Appendix B. The sharp-interface limit consists of asymptotic expansions in the bulk phases (in the variables $x$ and $y$), asymptotic expansions in the interface regions (in the variable $x$ and a new variable $z$ with characteristic length scale ${{\varepsilon }}$), and the matching of these asymptotic expansions.

6.2. The upscaled sharp-interface model

We will summarize the results of the matched asymptotic expansions presented in Appendix A. For this, we drop the subscript $0$ and the superscript $out$ for ease of notation. We call (6.1)–(6.19) the upscaled sharp-interface model.

The macroscopic equations for the unknowns $Q_f, p$ and $c$ are given by (B5), (B6) and (B9), that is

(6.1)$$\begin{gather} \partial_x Q_f = 0, \end{gather}$$

(6.2)$$\begin{gather}Q_f ={-}K_f \partial_x p, \end{gather}$$

(6.3)$$\begin{gather}\frac{{{\rm d}}}{{{\rm d}}t} ( \phi_{c,{total}} c ) + \partial_x ( ({-}K_c \partial_x p) c ) = \frac{1}{\overline{Pe_c}} \partial_x ( \phi_{c,{total}} \partial_x c) + \overline{Da} R_{interface}. \end{gather}$$

The coefficients of the upscaled equations depend on the distribution of the phases in the thin strip. In contrast to the upscaled phase-field model (5.12)–(5.27) the sharp-interface limit does not depend on the phase-field variables ${{\boldsymbol {\varPhi }}}$. Instead, the three disjoint domains $\varOmega _1(t), \varOmega _2(t)$ and $\varOmega _3(t)$ are used to locate the phases. The interface between $\varOmega _i$ and $\varOmega _j$ is denoted by $\varGamma _{ij}$. We introduce the notation $\varOmega _i|_x = \{{y\in [-\ell _\varOmega /2, \ell _\varOmega /2] : (x,y)\in \varOmega _i(t)}\}$, and write $N(\varGamma _{13})$ for the number of $\varGamma _{13}$ interfaces at a given $x$. With (B10), (B7), (B11) and (B53) we can calculate the coefficients of (6.1)–(6.3) as

(6.4)$$\begin{gather} \phi_{c,{total}}(x) = \text{vol}(\varOmega_1|_x), \end{gather}$$

(6.5)$$\begin{gather}K_f(t,x) = \int_{\varOmega_1|_x \cup \varOmega_2|_x} w\, {{\rm d}}y, \end{gather}$$

(6.6)$$\begin{gather}K_c(t,x) = \int_{\varOmega_1|_x} w\, {{\rm d}}y, \end{gather}$$

(6.7)$$\begin{gather}R_{interface} ={-} N(\varGamma_{13}) r(c). \end{gather}$$

We describe the evolution of the phases with the interface velocity $\nu$. This velocity in the $y$ direction is given by (B45), (B 39), (B40), summarized as

(6.8)$$\begin{gather} \nu ={\pm} \overline{Da} r(c) \text{ on } \varGamma_{13}, \text{ with } \varOmega_1 \text{ in the} \pm y \text{ direction}, \end{gather}$$

(6.9)$$\begin{gather}\nu = 0 \text{ on } \varGamma_{23}, \end{gather}$$

(6.10)$$\begin{gather}\nu ={-}(\partial_x s) {{\boldsymbol{v}}}^{(1)}_0 + {{\boldsymbol{v}}}^{(2)}_1 \text{ on } \varGamma_{12}. \end{gather}$$

For the flow profile we solve at each $x$ and $t$ a cell problem for the unknown $w$. Summarizing (B2), (B3), (B48), (B50) and the boundary condition (5.25), the unknown $w$ is given by the second-order differential equation

(6.11)$$\begin{gather} - \partial_y (\gamma_1 \partial_y w) = 1 \text{ in } \varOmega_1|_x, \end{gather}$$

(6.12)$$\begin{gather}- \partial_y (\gamma_2 \partial_y w) = 1 \text{ in } \varOmega_2|_x, \end{gather}$$

(6.13)$$\begin{gather}\rho_3 d_0 w - \partial_y (\gamma_3 \partial_y w) = 0 \text{ in } \varOmega_3|_x, \end{gather}$$

(6.14)$$\begin{gather}{[}\![ {w} ]\!] = 0 \text{ at } \varGamma_{12}, \varGamma_{13}, \varGamma_{23}, \end{gather}$$

(6.15)$$\begin{gather}{[}\![ {\gamma \partial_y w} ]\!] = 0 \text{ at } \varGamma_{12}, \varGamma_{13}, \varGamma_{23}, \end{gather}$$

(6.16)$$\begin{gather}w = 0 \text{ at } y ={\pm} \ell_\varOmega/2. \end{gather}$$

For the transport of the fluid–fluid interface $\varGamma _{12}$ in (6.10) we need the flow velocities ${{\boldsymbol {v}}}^{(1)}_0$ and ${{\boldsymbol {v}}}^{(2)}_1$. We then get the horizontal flow velocity ${{\boldsymbol {v}}}^{(1)}_0$ from (5.26), that is

(6.17)\begin{equation} {{\boldsymbol{v}}}^{(1)}_0 ={-} w \partial_x p_0. \end{equation}

For the vertical flow velocity ${{\boldsymbol {v}}}^{(2)}_1$ one has to solve (B4), (B29) and (B31), summarized

(6.18)$$\begin{gather} \partial_x ({{\boldsymbol{v}}}^{(1)}_{0}) + \partial_y ({{\boldsymbol{v}}}^{(2)}_{1}) = 0 \text{ in } \varOmega_1 \cup \varGamma_{12} \cup \varOmega_2, \end{gather}$$

(6.19)$$\begin{gather}-(\partial_x s) {{\boldsymbol{v}}}^{(1)}_0 + {{\boldsymbol{v}}}^{(2)}_{1} = 0 \text{ on } \varGamma_{13} \text{ and } \varGamma_{23}. \end{gather}$$

6.3. Upscaled sharp-interface model in a simplified geometry with symmetry

The upscaled sharp-interface model (6.1)–(6.19) uses no assumption on how the phases are distributed. When these are appearing in a fixed order, the model simplifies. In this case, there is no need to consider a general subdomain $\varOmega _i$ for the phase $i$, it is sufficient to know the width of the phase $i$ layer in the $y$ direction. These widths become unknowns of the model.

We assume here the following simplified geometry. The solid phase (in $\varOmega _3$) is covered by a film of fluid $1$ (occupying $\varOmega _1$). The second fluid (in $\varOmega _2$) is located in the middle of the thin strip. For simplicity, we assume symmetry around the $x$-axis. An illustration of the geometry is given in figure 2.

Figure 2. Symmetric geometry of two fluid phases in a thin strip.

With functions $d_1(t,x) > 0, d_2(t,x) > 0$, representing the width in the $y$ direction of the fluid phase 1, respectively 2, we can describe this situation by defining

(6.20)\begin{gather} \varOmega_2(t) = \{{(x,y) : -d_2(t,x) < y < d_2(t,x) }\}, \end{gather}

(6.21)\begin{gather} \varOmega_1(t) =\{{(x,y) : -d_1(t,x)-d_2(t,x) < y <{-}d_2(t,x)}\} \nonumber\\ \qquad\qquad\qquad \cup \{{(x,y) : d_2(t,x) < y < d_1(t,x)+d_2(t,x) }\} ,\end{gather}

(6.22)\begin{gather} \varOmega_3(t) =\{{(x,y) : -\ell_\varOmega/2 < y <{-}d_1(t,x)-d_2(t,x)}\} \nonumber\\ \qquad\qquad \qquad \cup \{{(x,y) : d_1(t,x)+d_2(t,x)< y < \ell_\varOmega/2 }\}. \end{gather}

In this geometry the solution $w$ to the cell problem (6.11)–(6.16) depends only on the variables $d_1$ and $d_2$, and on the choice of $\ell _\varOmega$. With a lengthy calculation we find that the terms depending on $\ell _\varOmega$ decay exponentially fast for big $\ell _\varOmega$, and we drop them in the following. The remaining terms lead to

(6.23)$$\begin{gather} K_f = \frac{2}{\gamma_1}\left(\frac{(d_1+d_2)^{3}}{3} + \left(\frac{\gamma_1}{\gamma_2} - 1\right)\frac{d_2^{3}}{3} + L_{slip}(d_1+d_2)^{2}\right), \end{gather}$$

(6.24)$$\begin{gather}K_c = \frac{2}{\gamma_1}\left(\frac{d_1^{3}}{3} + \frac{d_1^{2} d_2}{2} + L_{slip} d_1(d_1+d_2)\right), \end{gather}$$

with the slip length $L_{slip}$ given by

(6.25)\begin{equation} L_{slip} = \frac{\gamma_1}{\sqrt{\rho_3 d_0 \gamma_3}}. \end{equation}

We can relate $\partial _t d_1$ and $\partial _t d_2$ with the interface velocities (6.8)–(6.10). Considering the fluid–solid interface $\varGamma _{13}$ we get with (6.8)

(6.26)\begin{equation} \partial_t ( d_1 + d_2 ) = \nu ={-} \overline{Da} r(c), \end{equation}

while for the fluid–fluid interface $\varGamma _{12}$ we calculate with (6.10), (6.18) and (6.19)

(6.27)\begin{align} \partial_t d_2 &= \nu ={-} (\partial_x d_2) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_2) + {{\boldsymbol{v}}}^{(2)}_1(t,x,d_2) \nonumber\\ &={-} (\partial_x d_2) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_2) + {{\boldsymbol{v}}}^{(2)}_1(t,x,d_2) \nonumber\\ &\quad + (\partial_x (d_1 + d_2) ) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_1+d_2) - {{\boldsymbol{v}}}^{(2)}_1(t,x,d_1+d_2)\nonumber\\ &= (\partial_x (d_2+d_1)) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_1+d_2) - (\partial_x d_2) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_2) -\int_{d_2}^{d_2+d_1} \partial_y {{\boldsymbol{v}}}^{(2)}_1(t,x,y)\, {{\rm d}}y \nonumber\\ &= (\partial_x (d_2+d_1)) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_1+d_2) - (\partial_x d_2) {{\boldsymbol{v}}}^{(1)}_0(t,x,d_2) +\int_{d_2}^{d_2+d_1} \partial_x {{\boldsymbol{v}}}^{(1)}_0(t,x,y)\, {{\rm d}}y \nonumber\\ &= \partial_x \left( \int_{d_2}^{d_2+d_1} {{\boldsymbol{v}}}^{(1)}_0(t,x,y)\, {{\rm d}}y \right) . \end{align}

The integral equals the total fluid flux in the $x$ direction in the upper half of $\varOmega _1$. We use (6.17), (6.6) and the symmetry of $w$ around $y=0$ to further calculate

(6.28)\begin{equation} \partial_t d_2 = \partial_x \left( \int_{d_2}^{d_2+d_1} {{\boldsymbol{v}}}^{(1)}_0\, {{\rm d}}y \right) ={-}\partial_x \left( (\partial_x p) \int_{d_2}^{d_2+d_1} w\, {{\rm d}}y \right) ={-}\frac{1}{2} \partial_x\left(K_c \partial_x p\right). \end{equation}

We can now summarize (6.1), (6.2), (6.3), (6.26) and (6.26) as an upscaled model for the unknowns $d_1, d_2, p, Q_f$ and $c$

(6.29)$$\begin{gather} \partial_t d_1 + \partial_t d_2 ={-} \overline{Da} r(c(t,x)), \end{gather}$$

(6.30)$$\begin{gather}\partial_t d_2 ={-} \frac{1}{2} \partial_x (K_c(d_1,d_2) \partial_x p), \end{gather}$$

(6.31)$$\begin{gather}Q_f ={-}K_f(d_1,d_2) \partial_x p, \end{gather}$$

(6.32)$$\begin{gather}\partial_x Q_f = 0, \end{gather}$$

(6.33)$$\begin{gather}\frac{{{\rm d}}}{{{\rm d}}t} ( 2 d_1 c ) + \partial_x ( ({-}K_c(d_1,d_2) \partial_x p) c ) = \frac{1}{Pe_c} \partial_x ( 2 d_1 \partial_x c) - 2 \overline{Da} r(c) . \end{gather}$$

Remark 6.1 We can rewrite (6.30) and (6.31) to highlight the hyperbolicity of the model. As discussed in Remark 5.1 one assumption for the upscaling is that there is no occurrence of triple points. Therefore, we assume $d_1 > 0$ and $d_2 > 0$ and deduce $K_f > 0, K_c > 0$. We can now calculate

(6.34)\begin{equation} \partial_t d_2 = \frac{1}{2} Q_f \partial_x \left(\frac{K_c(d_1,d_2)}{K_f(d_1,d_2)} \right). \end{equation}

The unknown $d_2$ gets transported with flux $Q_f K_c / K_f$ and can show hyperbolic behaviour, such as the formation of discontinuities.

6.4. Asymptotic consistency

In § 5 we have investigated the limit process $\beta \to 0$, while in § 6 we examined $\bar {{\varepsilon }} \to 0$. A common question is under which circumstances there is asymptotic consistency, i.e. these two limit processes commute. In figure 3 all limit processes are shown in a commutative diagram.

Figure 3. Models obtained by upscaling ($\beta \to 0$) and the sharp-interface limit ($\bar {{\varepsilon }} \to 0$).

We investigate asymptotic consistency with non-dimensional numbers chosen as in (5.5)–(5.10) with $\overline {Re}, \overline {Ca}, \overline {Da}, \bar M, \overline {Pe_c}$ constant and independent of $\bar {{\varepsilon }}$ and $\beta$. The non-dimensional $\delta$ is chosen as $\delta = \bar {{\varepsilon }}$.

When starting with the fully resolved diffuse-interface model (4.3)–(4.10) the limit $\bar {{\varepsilon }} \to 0$ results in a sharp-interface model as described in § 2. For details on this sharp-interface limit, see Rohde & von Wolff (2021).

When we assume the geometry of § 6.3 we can proceed to upscale the fully resolved sharp-interface model after introducing $d_1$ and $d_2$. While the process is tedious, the main ideas are analogous to the calculations in Sharmin et al. (2020). In particular, the asymptotic expansion of interface conditions, normal vectors and curvature have to be handled with care, as the coordinates ${{\boldsymbol {x}}} = (x,\beta y)$ depend on $\beta$. For sake of brevity we skip this calculation here.

With the geometry of § 6.3 we find asymptotic consistency, that is the limit processes $\beta \to 0$ and $\bar {{\varepsilon }} \to 0$ commute. The result of the upscaling of the fully resolved sharp-interface model is exactly given by (6.29)–(6.33).

Remark 6.2 In more general geometries, asymptotic consistency does not necessary hold. This is due to the following observation. When upscaling a fully resolved diffuse-interface model, the parameter $\delta$ is constant and independent of $\beta$. This leads to $\tilde \phi _f > 0$ and $\tilde \phi _c > 0$ everywhere. Because of this, we obtain upscaled equations for $p$ and $c$ without further assumptions on the geometry. The upscaled variables $p$ and $c$ do not depend on $y$, even if the geometry consists of two parallel channels separated by a solid region with ${{\boldsymbol {\varPhi }}} \approx {{\boldsymbol {e}}}_3$. On the other hand, when upscaling the fully resolved diffuse-interface model, the $\delta$-modifications have already vanished, as $\delta = \bar {{\varepsilon }}$. In this case, it is possible to have a different pressure $p$ in each channel, that is in each connected part of $ {\varOmega _1|_x {\bar \cup} \varOmega _2|_x}$. Also, it is possible to have a different ion concentration $c$ in each connected part of $\varOmega _1|_x$.

We conclude that we have asymptotic consistency under the condition that there is only one flow channel, i.e. ${\varOmega _1|_x {\bar \cup} \varOmega _2|_x}$ is connected for every $x$, and that the first fluid phase is connected, i.e. $\varOmega _1|_x$ is connected for every $x$. It is also possible to consider the symmetric case as in § 6.3 and have two symmetric connected parts of fluid one.

7. Numerical investigation

We will now compare the upscaled $\delta$-$2f1s$-model (5.12)–(5.27) with the fully resolved $\delta$-$2f1s$-model (4.3)–(4.10). Remark 6.1 suggests that shock fronts can form in the upscaled model. Note that in this case the assumptions for the upscaling in § 5 are no longer valid, and we expect different behaviours from the two models.

For the fully resolved $\delta$-$2f1s$-model (4.3)–(4.10) we use a monolithic finite-element implementation provided by the DUNE-Phasefield module (von Wolff 2021). We employ Taylor–Hood elements for the flow variables velocity and pressure, and first-order Lagrange elements for the ion concentration and the phase-field parameters. The implementation is based on DUNE-PDELab (Bastian, Heimann & Marnach 2010) using ALU-Grid routines for adaptive grid generation (Alkämper et al. 2016).

7.1. Numerical scheme for the upscaled $\delta$-$2f1s$-model

The upscaled $\delta$-$2f1s$-model consists of multiple coupled problems. The upscaled (5.12)–(5.14) for the unknowns $Q_f, p$ and $c$ have parameters (5.15)–(5.18) that depend on the distribution of phases in the $y$-direction. This distribution is described by the fully coupled two-dimensional problem (5.19)–(5.22) for the Cahn–Hilliard variables $\phi _1, \phi _2, \phi _3, \mu _1, \mu _2, \mu _3$. Furthermore, the flow profile has to be calculated by the cell problem (5.24) and (5.25).

For simplicity, we present the numerical scheme for equidistant time steps $t_n = n {\rm \Delta} t$ and equidistant discretization in $x$ by $x_k = k {\rm \Delta} x$. Let also $x_{k+1/2} = (x_k + x_{k+1})/2$. For each $t_n, x_k$ we discretize the one-dimensional unknown $\phi ^{n}_{1,k}(y) = \phi _1(t_n,x_k,y)$ with linear Lagrange elements, and analogously for $\phi ^{n}_{2,k}, \phi ^{n}_{3,k}, \mu ^{n}_{1,k}, \mu ^{n}_{2,k}, \mu ^{n}_{3,k}, {{\boldsymbol {v}}}^{(1),n}_{0,k}, {{\boldsymbol {v}}}^{(2),n}_{1,k}, w^{n}_k$. Again, we also use this notation for other variables such as $\tilde \phi ^{n}_{f,k}$.

We discretize the macroscopic unknown $c(t,x)$ with a finite volume scheme, that is $c^{n}(x) = c(t_n,x)$ is piecewise constant with $c(t_n,x) = c^{n}_k$ for $x\in (x_{k-1/2},x_{k+1/2})$. The pressure $p^{n}(x) = p(t_n,x)$ is discretized using linear Lagrange elements with nodes $x_{k+1/2}$. Therefore, $\partial _x p$ is constant on each finite volume cell $(x_{k-1/2},x_{k+1/2})$.

Given ${{\boldsymbol {\varPhi }}}^{n}_k, c^{n}_k$ for all $x_k$ at time $t_n$, we now calculate the next time step using the following algorithm.

1. (i) For each $x_k$ use (5.24) and (5.25) to solve for $w^{n}_k(y)$. Here, we use ${{\boldsymbol {\varPhi }}} = {{\boldsymbol {\varPhi }}}^{n}_k$ and the finite element method to discretize the equation. The equations for different $x_k$ are independent and can be solved in parallel.

2. (ii) For each $x_k$ calculate $K^{n}_{f,k}$ and $K^{n}_{c,k}$ by

   (7.1a,b)\begin{equation} K^{n}_{f,k} = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \tilde \phi^{n}_{f,k} w^{n}_k\, {{\rm d}}y, \quad K^{n}_{c,k} = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \tilde \phi^{n}_{c,k} w^{n}_k\, {{\rm d}}y. \end{equation}

3. (iii) Solve for $p^{n}(x)$ using the finite element method with

   (7.2)\begin{equation} \partial_x( - K^{n}_f \partial_x p^{n}) = 0. \end{equation}
   Here, $K^{n}_f(x) = K^{n}_{f,k}$ for $x\in (x_{k-1/2},x_{k+1/2})$. As $K^{n}_f > 0$, the pressure $p$ is either a monotonically increasing or monotonically decreasing function, depending on the boundary conditions. We assume from here on $\partial _x p^{n} \leq 0$ and therefore fluid flow in the positive $x$ direction. In case $\partial _x p^{n} \geq 0$ the upwind schemes in steps (v) and (vii) have to be modified.

4. (iv) For each $x_k$ calculate ${{\boldsymbol {v}}}^{(1),n}_{0,k}(y) = - w^{n}_k(y) \partial _x p^{n}(x_k)$.

5. (v) Next, for each $x_k$ we solve for ${{\boldsymbol {v}}}^{(2),n}_{1,k}$ and the Cahn–Hilliard variables $\phi ^{n+1}_{2,k}, \phi ^{n+1}_{3,k}, \mu ^{n+1}_{1,k}, \mu ^{n+1}_{2,k}, \mu ^{n+1}_{3,k}$. For ${{\boldsymbol {v}}}^{(2),n}_{1,k}$ we use (5.27) the with an explicit upwind scheme for the $x$-derivative, i.e.

   (7.3)\begin{equation} \partial_y (\tilde \phi^{n+1}_{f,k} {{\boldsymbol{v}}}^{(2),n}_{1,k}) ={-} \frac{\tilde \phi^{n}_{f,k} {{\boldsymbol{v}}}^{(1),n}_{0,k} - \tilde \phi^{n}_{f,k-1} {{\boldsymbol{v}}}^{(1),n}_{0,k-1} }{{\rm \Delta} x}. \end{equation}
   This equation is coupled with the Cahn–Hilliard cell problems (5.19)–(5.22). We again use an explicit upwinding scheme for the $x$-derivative, that is
   (7.4)\begin{gather} \frac{\phi^{n+1}_{1,k}-\phi^{n}_{1,k}}{{\rm \Delta} t} + \frac{\phi^{n}_{1,k} {{\boldsymbol{v}}}^{(1),n}_{0,k} - \phi^{n}_{1,k-1} {{\boldsymbol{v}}}^{(1),n}_{0,k-1}}{{\rm \Delta} x} + \partial_y (\phi^{n+1}_{1,k} {{\boldsymbol{v}}}^{(2),n}_{1,k}) - \frac{\bar {{\varepsilon}} \bar M}{\varSigma_1} \partial_y^{2} \mu^{n+1}_{1,k} \nonumber\\ \hspace{-4pc}={-} \frac{\overline{Da}}{\bar {{\varepsilon}} } q({{\boldsymbol{\varPhi}}}^{n+1}_k)\left(r(c^{n}(x_k)) + \tilde \alpha \mu^{n+1}_{1,k} - \tilde \alpha \mu^{n+1}_{3,k}\right), \end{gather}
   (7.5)$$\begin{gather} \frac{\phi^{n+1}_{2,k}-\phi^{n}_{2,k}}{{\rm \Delta} t} + \frac{\phi^{n}_{2,k} {{\boldsymbol{v}}}^{(1),n}_{0,k} - \phi^{n}_{2,k-1} {{\boldsymbol{v}}}^{(1),n}_{0,k-1}}{{\rm \Delta} x} + \partial_y (\phi^{n+1}_{2,k} {{\boldsymbol{v}}}^{(2),n}_{1,k}) - \frac{\bar {{\varepsilon}} \bar M}{\varSigma_1} \partial_y^{2} \mu^{n+1}_{2,k} = 0, \end{gather}$$
   (7.6)$$\begin{gather}\phi^{n+1}_{3,k} = 1- \phi^{n+1}_{1,k} -\phi^{n+1}_{2,k}, \end{gather}$$
   (7.7)$$\begin{gather}\mu^{n+1}_{1,k} = \frac{\partial_{\phi_1} W({{\boldsymbol{\varPhi}}}^{n+1}_k)}{\bar {{\varepsilon}}} - \bar {{\varepsilon}} \varSigma_i \partial_y^{2} \phi^{n+1}_{1,k}, \end{gather}$$
   (7.8)$$\begin{gather}\mu^{n+1}_{2,k} = \frac{\partial_{\phi_2} W({{\boldsymbol{\varPhi}}}^{n+1}_k)}{\bar {{\varepsilon}}} - \bar {{\varepsilon}} \varSigma_i \partial_y^{2} \phi^{n+1}_{2,k}, \end{gather}$$
   (7.9)$$\begin{gather}\mu^{n+1}_{3,k} ={-} \mu^{n+1}_{1,k} - \mu^{n+1}_{2,k}. \end{gather}$$
   Note that we do not use (5.21) and (5.22) for $\phi ^{n+1}_{3,k}$ and $\mu ^{n+1}_{3,k}$. Instead, we use that by construction $\phi _1 + \phi _2 + \phi _3 = 1$ and $\mu _1 + \mu _2 + \mu _3 = 0$, see (Rohde & von Wolff 2021) for details.

   We use the finite element method to discretize (7.3)–(7.9) and Newtons method to solve the resulting nonlinear system. This step has by far the highest computational cost. With the explicit upwinding scheme for the $x$ derivatives, the cell problems for each $k$ fully decouple and can be solved in parallel. This leads to a significant speed up in comparison with discretizing the Cahn–Hilliard evolution (5.19)–(5.22) naively as a two-dimensional problem.

6. (vi) Calculate $\tilde \phi ^{n+1}_{c,{total},k}$ and $R^{n+1}_{{total},k}$ as

   (7.10)$$\begin{gather} \tilde \phi^{n+1}_{c,{total},k} = \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} \phi^{n+1}_{c,k}\, {{\rm d}}y, \end{gather}$$
   (7.11)$$\begin{gather}R^{n+1}_{{total},k} ={-} \int_{-\ell_\varOmega/2}^{\ell_\varOmega/2} q({{\boldsymbol{\varPhi}}}^{n+1}_k)(r(c^{n}(x_k)) + \tilde \alpha \mu^{n+1}_{1,k} - \tilde \alpha \mu^{n+1}_{3,k})\,{{\rm d}}y. \end{gather}$$
   We also set $\tilde \phi ^{n+1}_{c,{total},k+1/2} = (\tilde \phi ^{n+1}_{c,{total},k} +\tilde \phi ^{n+1}_{c,{total},k+1})/2$.

7. (vii) Finally we solve for $c$ using (5.14) discretized by the finite volume method. We use an implicit upwinding scheme for the transport in the $x$-direction

   (7.12)\begin{align} &\frac{\tilde \phi^{n+1}_{c,{total},k} c^{n+1}_k - \tilde \phi^{n}_{c,{total},k} c^{n}_k}{{\rm \Delta} t} - \frac{ K^{n}_{c,k} \partial_x p^{n}(x_k) c^{n+1}_k - K^{n}_{c,k-1} \partial_x p^{n}(x_{k-1}) c^{n+1}_{k-1}}{{\rm \Delta} x} \nonumber\\ &\quad = \frac{1}{\overline{Pe_c}} \frac{1}{{\rm \Delta} x} \left( \tilde \phi^{n+1}_{c,{total},k+1/2} \frac{c^{n+1}_{k+1} - c^{n+1}_{k}}{{\rm \Delta} x} - \phi^{n+1}_{c,{total},k-1/2} \frac{c^{n+1}_k - c^{n+1}_{k-1}}{{\rm \Delta} x}\right) + \frac{\overline{Da}}{\bar {{\varepsilon}} } R^{n+1}_{{total},k}. \end{align}

7.2. Comparison: formation of an $N$-wave

As our first numerical example we choose a geometry as described in § 6.3, with the computational domain $(x,y) \in [0,1]\times [-1,0]$. For $x = 0$ and $x = 1$ we choose periodic boundary conditions for all variables except the pressure $p$. For $y = -1$ we use the trivially upscaled versions of the boundary conditions 5.1-5.4 and for $y=0$ we choose boundary conditions according to the symmetry assumption.

We will compare the non-dimensional $\delta$-$2f1s$-model with the upscaled $\delta$-$2f1s$-model (5.12)–(5.27). For simplicity, we choose $\gamma _1 = \gamma _2$ and $d_0$ sufficiently big such that $L_{slip} \approx 0$. We choose the phase-field parameter $\bar {{\varepsilon }} = 0.03$ and $\delta = \bar {{\varepsilon }}$ as in § 6.

We want to focus on the hyperbolic behaviour of $d_2$ as described in Remark 6.1. Therefore, we choose $c$ in the initial conditions such that $r(c) = 0$. This leads to no precipitation or dissolution in the model, and the fluid–solid interface does not change over time. We choose

(7.13a,b)\begin{equation} d_1 + d_2 \equiv 0.7\quad \text{and}\quad d_1(x) = 0.4 + 0.15\sin(2{\rm \pi} x). \end{equation}

This corresponds to a plane fluid–solid interface and a sine-shaped fluid–fluid interface. An image of these initial conditions is given in figure 4.

Figure 4. Evolution of the upscaled $\delta$-$2f1s$-model on the domain $[0,1]\times [-1,0]$. Shown in red is fluid phase one, with fluid phase two above and solid phase below. From left to right: initial data, $t=0.15, t=0.3$ and $t=0.45$.

By applying a pressure difference as Dirichlet boundary condition at $x=0$ and $x=1$, the two fluid phases will move in the positive $x$-direction. The fluid velocity ${{\boldsymbol {v}}}^{(1)}_0$ is higher in the centre of the channel. As shown in figure 4, this will lead to a steeper fluid–fluid interface over time. At a time $t^{\ast }>0$ the upscaled $\delta$-$2fs$-model has a fluid–fluid interface that is perpendicular to the thin strip. As discussed in Remark 5.1, the assumptions for the upscaling in § 5 are no longer valid. For times $t>t^{\ast }$ the fluid–fluid interface will roll over, leading to multiple layers of fluid phase 1 at the same $x$ value. In § 8 we will present some ideas of how to handle such cases in future works.

We can compare this behaviour with the non-dimensional $\delta$-$2f1s$-model in a thin strip for different values of $\beta$. As shown in figure 5, for times $t < t^{\ast }$ there is a good agreement between the non-dimensional $\delta$-$2f1s$-model with small values of $\beta$ and the upscaled $\delta$-$2f1s$-model.

Figure 5. Fluid–fluid interface locations for the non-dimensional $\delta$-$2f1s$-model with varying $\beta$, and for the upscaled $\delta$-$2f1s$-model. The interface is located through the condition $\phi _1 = \phi _2$. Left: $t=0.3$, right: $t=0.44$.

In contrast to the upscaled $\delta$-$2f1s$-model, the non-dimensional $\delta$-$2f1s$-model does not evolve to a fluid–fluid interface perpendicular to the thin strip, as shown in figure 5. Instead, when reaching a steep fluid–fluid interface there are regions of high curvature at the beginning and end of the steep passage. In these regions of high curvature the surface tension leads to a pressure difference between the fluid phases, counteracting the interface getting steeper. For smaller $\beta$ the fluid–fluid interface allows for a steeper passage in $(x,y)$ coordinates, as this effect depends on the curvature in the ${{\boldsymbol {x}}}$ coordinates, which are not scaled with $\beta$.

7.3. Comparison: precipitation

In the second numerical example we study precipitation in the thin strip. We use the same domain and boundary conditions as in the previous example. Again, we choose $\gamma _1 = \gamma _2$, and a $d_0$ large enough so that $L_{slip}\approx 0$. We further choose $\bar {{\varepsilon }} = 0.03$ and $\delta = \bar {{\varepsilon }}$. We use a simple, linear reaction rate $r(c) = c - 0.5$ and choose the ion concentration to be in equilibrium initially, that is $c = 0.5$ everywhere. With $d_1(x) = 0.4$ and $d_2(x) = 0.3$ in the initial conditions correspond to the phases being layered in the thin strip, without depending on $x$. To induce precipitation we add a source term $s(x)$ to the ion conservation equation (4.5), it now reads

(7.14)\begin{equation} \partial_{t} (\tilde\phi_c c) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\phi_c {{\boldsymbol{v}}}c) + \frac{Cn}{\beta Pe_{CH}} \boldsymbol{\nabla} \boldsymbol{\cdot}({{\boldsymbol{J}}}_1 c) = \frac{1}{Pe_c} \boldsymbol{\nabla} \boldsymbol{\cdot} ( \tilde\phi_c \boldsymbol{\nabla} c) + Da R + \tilde \phi_c s(x). \end{equation}

The source terms upscales trivially at $O(\beta ^{0})$, and the upscaled ion conservation equation (5.14) is now given by

(7.15)\begin{equation} \frac{{{\rm d}}}{{{\rm d}}t} ( \tilde \phi_{c,{total}} c ) + \partial_x ( ({-}K_c \partial_x p) c ) = \frac{1}{\overline{Pe_c}} \partial_x ( \tilde \phi_{c,{total}} \partial_x c) + \frac{\overline{Da}}{\bar {{\varepsilon}} } R_{total} + \tilde \phi_{c,{total}} s(x). \end{equation}

We choose the ion source to be located between $x = 0.1$ and $x=0.3$, in detail

(7.16)\begin{equation} s(x) = \max(0, 62.5 (x-0.1)(0.3-x)). \end{equation}

Figure 6 shows a comparison between the non-dimensional $\delta$-$2f1s$-model with different values of $\beta$, and the upscaled $\delta$-$2f1s$-model. There is a good agreement between the full model with small values of $\beta$ and the upscaled model. For large values of $\beta$ there is less precipitation in the thin strip. This is due to the ion concentration $c$ not being constant in the $y$-direction. The source term $\tilde \phi _c s(x)$ generates ions everywhere in the first fluid phase, but precipitation removes ions from the first fluid phase only at the fluid–solid interface. This leads to an oversaturation $c>0.5$ further away from the fluid–solid interface. For smaller values of $\beta$ the diffusion in the $y$-direction results in more ions precipitating and therefore a smaller oversaturation of ions in the fluid phase.

Figure 6. Interface locations at time $t=2.4$ for the non-dimensional $\delta$-$2f1s$-model with varying $\beta$, and for the upscaled $\delta$-$2f1s$-model. The fluid–fluid interface can be seen in the upper half and is located by the condition $\phi _1 = \phi _2$. The fluid–solid interface in the lower half is located by $\phi _1 = \phi _3$. The initial location of the fluid–solid interface is marked by a black line.

Figure 6 also shows the influence of a non-constant width of the thin strip on the flow inside the thin strip. The fluid–fluid interfaces are pushed towards the centre of the thin strip, where flow velocities are higher.

8. Conclusion

We have upscaled a phase-field model for the incompressible flow of two immiscible fluids in the geometry of a thin strip. With the assumption of slow variation in the direction along the thin strip, we find a two-scale model as a result of the upscaling. The two-scale model consists of macroscopic equations for the total flux $Q_f$, the pressure $p$ and the ion concentration $c$. Those equations are coupled to microscopic equations for the flow and for the geometry, represented by the phase-field variables. By discretizing the upscaled equations with a finite volume scheme, we obtain microscopic cell problems that are fully decoupled in each time step and can therefore be solved in parallel.

We have also investigated the sharp-interface limit of the upscaled phase-field model and found under additional assumptions on the geometry a fully upscaled model. This model only consists of macroscopic equations for total flux $Q_f$, the pressure $p$, the ion concentration $c$ and the widths of each fluid phase, $d_1$ and $d_2$. Further analysis shows that the upscaling and the sharp-interface limit commute.

Further research is needed to generalize this work in a multitude of directions. The most obvious next step is to consider three-dimensional geometries, such as a thin tube. In such a case, the curvature of the fluid–fluid interface will be bigger by a factor of $\beta ^{-2}$ compared with the two-dimensional case. With non-dimensional numbers scaling as in this work, surface tension effects will enter the leading-order equations of the momentum equation. As a result, one has to introduce a faster time scale to resolve relaxation in the cross-section towards an equilibrium state.

A second possible direction for future research is the consideration of a finite number of points where the assumption of slow variation along the thin strip is not fulfilled. This includes three-phase contact points, as well as the $N$ waves shown in the numerical investigation. Such cases can not be upscaled with the current assumptions, and the numerical investigation shows no agreement between the fully resolved and the upscaled model. In future models it might be possible to describe regions with fast variation along the thin strip as fully resolved, and couple these regions on each side with the upscaled phase-field model.

Lastly, by increasing the complexity of the underlying phase-field model, the upscaling procedure will become increasingly difficult. Interesting effects seen in experiments include the dissolved ions being driven by electrostatic fields, as well as small mineral particles getting transported by the fluid flow.