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 (Reference Rohde and von Wolff2021), 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 (Reference Rohde and von Wolff2021) 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.

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. Reference Sharmin, Bringedal and Pop2020; Lunowa et al. Reference Lunowa, Bringedal and Pop2021). Note that this is a major difference to the three-dimensional case (see e.g. Mikelić Reference 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.

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.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 (Reference Pego and Penrose1989).

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 (Reference Rohde and von Wolff2021).

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. (Reference Sharmin, Bringedal and Pop2020). 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).

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 Reference Rohde and von Wolff2021) 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.