2.1. Starting assumptions and system configuration

The incompressible, inertial and quasi-steady flow under consideration is that of two non-Newtonian fluid phases, $\beta$ and $\gamma$, through a rigid and homogeneous porous medium (for which the solid skeleton is denoted as the $\sigma$ phase). By quasi-steady it is meant that the shapes and positions of the fluid–fluid interfaces within the pores are not stationary, although momentum transfer in the bulk of each phase is considered steady. Furthermore, the derivations that follow are applicable in situations in which both fluid phases circulate in either a continuous or discontinuous manner. Under these conditions, the mass and momentum balance equations at the pore scale are written as (in the following, $\alpha =\beta,\gamma$ is used to represent either one of the fluid phases)

(2.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_\alpha=0, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(2.1b)$$\begin{gather}\rho_\alpha\boldsymbol{v}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v}_\alpha=\rho_\alpha\boldsymbol{b}_\alpha + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{p_\alpha} \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

where ${\boldsymbol{\mathsf{T}}}_{p_\alpha }$ is the total stress tensor in the $\alpha$ phase, defined by

(2.1c)\begin{equation} {\boldsymbol{\mathsf{T}}}_{p_\alpha}={-}{\boldsymbol{\mathsf{I}}}p_\alpha+\mu(\varGamma_\alpha) ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+\boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ). \end{equation}

In these equations, $\boldsymbol {v}_\alpha$, $\rho _\alpha$ and $p_\alpha$ are, respectively, the velocity, density and pressure of the $\alpha$ phase, whereas ${\boldsymbol{\mathsf{I}}}$ is the identity tensor and $\boldsymbol {b}_\alpha$ represents the body force per unit mass, which is assumed to be constant within $\mathscr {V}_\alpha$. Moreover, in the viscous part of the total stress tensor, a generalised Newton's law of viscosity is used, which means that the apparent viscosity, $\mu (\varGamma _\alpha )$, is a function of the shear-rate modulus, $\varGamma _\alpha$, defined as

(2.1d)\begin{equation} \varGamma_\alpha =\sqrt{ \tfrac{1}{2} [ ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+ \boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ): ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+ \boldsymbol{\nabla} \boldsymbol{v}^T_\alpha )]}. \end{equation}

While writing (2.1b), it is assumed that, albeit convective acceleration may be significant in accordance with a Reynolds number not exceedingly small compared with unity, the frequency parameter in the $\alpha$ phase, defined as $\rho _\alpha \ell _{{ref} \alpha }^2/(\mu _{{ref} \alpha } t_{ref})$, remains small compared with 1 ($\ell _{{ref} \alpha }$, $\mu _{{ref} \alpha }$ and $t_{ref}$ being, respectively, a reference length and viscosity in the $\alpha$ phase and a reference time; these definitions are further discussed in § 6 once a particular rheological model is chosen).

2.2. Flow problem in a periodic unit cell

Following a classical assumption in upscaling, the analysis is restricted to systems in which there is a separation of length scales. This means that the following developments correspond to two-phase flow far away from a displacement front, if present, or within the entire system if the two phases coexist mixed as exemplified in figure 1. In the former case, it is assumed that the porous medium (of characteristic length $\mathscr {L}$) can be subdivided into regions (of characteristic length $L<\mathscr {L}$) that have a hierarchical nature. In these circumstances, an averaging domain, $\mathscr {V}$ (of measure $V$ and characteristic size $r_0\ll L$), representative of the pore-scale configuration and transport processes (at least for momentum as discussed below) can be assumed to exist. For this reason, $r_0$ must also be much larger than the characteristic pore-scale length, $\max (\ell _\beta,\ell _\gamma )$. This concept, which is introduced in the form of a periodic unit cell at the beginning in the homogenisation technique (Auriault, Boutin & Geindreau Reference Auriault, Boutin and Geindreau2009) or at the closure level in the volume-averaging method (Whitaker Reference Whitaker1999), is necessary for a local description. At this point, it is pertinent to mention that the periodicity assumption is a commodity more than a necessity in the upscaling approach used here. It allows the closure problems to be solved only in a unit cell, thus leading to local effective-medium coefficients. Furthermore, the upscaled models can be used in practice even in situations when the real porous medium is not periodic. In other words, periodicity is a convenient assumption for the closure problem definition and it has a minor impact at the macroscopic level. In the following, the averaging domain is denoted as $\mathscr {V}$ and the space occupied by each fluid phase is $\mathscr {V}_\alpha$ (of measure $V_\alpha$).

Figure 1. (a) Sketch of a porous medium saturated with two fluid phases, $\beta$ (in blue) and $\gamma$ (in white), flowing in the presence of a displacement front or as coexisting mixed phases including the scales, averaging domain and characteristic lengths. (b) Example of a periodic unit cell, of side length $\ell$, illustrating the position vectors locating the barycentres of each fluid phase ($\boldsymbol {x}_\beta$ and $\boldsymbol {x}_\gamma$), as well as an example of the position vectors that locate the fluid–fluid interface with respect to a fixed system of coordinates ($\boldsymbol {r}_{\beta \gamma }$) and with respect to the phase barycentres ($\boldsymbol {z}_{\beta \gamma }$ and $\boldsymbol {z}_{\gamma \beta }$).

In terms of the averaging domain, the superficial and intrinsic averages of a pore-scale quantity, $\psi _\alpha$, defined at any point within the $\alpha$ phase, located by $\boldsymbol {r}_{\alpha }$ with respect to a fixed system of coordinates, are respectively defined by the following integral operators:

(2.2a)$$\begin{gather} \langle\psi_\alpha\rangle_\alpha|_{\boldsymbol{x}_\alpha}= \frac{1}{V}\int_{\mathscr{V}_\alpha} \psi_\alpha|_{\boldsymbol{r}_\alpha}\, {\rm d} V, \end{gather}$$

(2.2b)$$\begin{gather}\langle \psi_\alpha \rangle^\alpha |_{\boldsymbol{x}_\alpha} = \varepsilon_\alpha^{{-}1} \langle \psi_\alpha \rangle_\alpha |_{\boldsymbol{x}_\alpha} = \frac{1}{V_\alpha}\int_{\mathscr{V}_\alpha} \psi_\alpha|_{\boldsymbol{r}_\alpha}\, {\rm d} V. \end{gather}$$

The two resulting averages are assigned at the barycentre of each phase, $\boldsymbol {x}_\alpha =\langle \boldsymbol {r}_\alpha \rangle ^\alpha$, of $\mathscr {V}_\alpha$. The volume fraction of the $\alpha$ phase is denoted by $\varepsilon _\alpha =V_\alpha /V=\varepsilon S_\alpha$, where $\varepsilon =(V_\beta +V_\gamma )/V$ is the porosity of the porous medium and $S_\alpha$ the $\alpha$ phase saturation. From the above definition, a pore-scale quantity can be expressed in terms of its intrinsic average and deviation by means of the following decomposition (Gray Reference Gray1975):

(2.3)\begin{equation} \psi_\alpha |_{\boldsymbol{r}_\alpha} = \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{r}_\alpha} + \tilde{\psi}_\alpha |_{\boldsymbol{r}_\alpha}.\end{equation}

Furthermore, the following Taylor series expansion can be used in order to evaluate $\langle \psi _\alpha \rangle ^\alpha$ at the phase barycentre:

(2.4)\begin{equation} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{r}_\alpha} =\langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha} +\boldsymbol{z}_\alpha\boldsymbol{\cdot} \boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha} +\tfrac{1}{2} \boldsymbol{z}_\alpha\boldsymbol{z}_\alpha: \boldsymbol{\nabla}\boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha}+\cdots ,\end{equation}

where $\boldsymbol {z}_\alpha =\boldsymbol {r}_\alpha -\boldsymbol {x}_\alpha$ locates the same point as $\boldsymbol {r}_\alpha$ but with respect to $\boldsymbol {x}_\alpha$ as illustrated in figure 1(b). Note that, by definition, $\langle \boldsymbol {z}_\alpha \rangle ^\alpha = \boldsymbol {0}$ and $\boldsymbol {z}_\alpha = \boldsymbol {O}(r_0)$. If this expansion is used for $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$, since this quantity can be assumed to be slowly varying (i.e. its spatial variations are characterised by $L$), it is reasonable to assume, on the basis of the length-scale constraint $r_0\ll L$, that

(2.5)\begin{equation} \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha|_{\boldsymbol{r}_\alpha}\simeq \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha|_{\boldsymbol{x}_\alpha}\equiv \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha, \end{equation}

which means that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ can be considered as constant at the scale $r_0$.

The pore-scale momentum balance equation is hence written in a periodic unit cell under the following form:

(2.6)\begin{equation} \rho_\alpha\boldsymbol{v}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}_\alpha={-}\boldsymbol{\nabla} \langle p_\alpha\rangle^\alpha +\rho_\alpha\boldsymbol{b}_\alpha + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha} \quad\text{in }\mathscr{V}_\alpha,\end{equation}

with

(2.7)\begin{equation} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha}={-}{\boldsymbol{\mathsf{I}}} \tilde{p}_\alpha+\mu(\varGamma_\alpha) ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+\boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ). \end{equation}

To complete the statement of the flow problem, it is assumed that no mass transport and no slip take place at the solid–fluid interfaces, $\mathscr {A}_{\alpha \sigma }$, which leads to the following boundary condition:

(2.8a)\begin{equation} \boldsymbol{v}_\alpha=\boldsymbol{0}, \quad \text{at }\mathscr{A}_{\alpha \sigma}. \end{equation}

Moreover, the same assumption is imposed at the fluid–fluid interface, $\mathscr {A}_{\beta \gamma }$. Hence, the mass balance boundary condition at this interface can be written as

(2.8b)\begin{equation} \boldsymbol{v}_\beta= \boldsymbol{v}_\gamma = \boldsymbol{w} \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{equation}

where $\boldsymbol {w}$ represents the speed of displacement of $\mathscr {A}_{\beta \gamma }$. In this work, no triple-phase contact is considered. Nevertheless, the pore-scale flow model statement considered here has been used by many authors (e.g. Auriault & Sanchez-Palencia Reference Auriault and Sanchez-Palencia1986; Whitaker Reference Whitaker1986, Reference Whitaker1994; Lasseux, Quintard & Whitaker Reference Lasseux, Quintard and Whitaker1996; Lasseux, Ahmadi & Arani Reference Lasseux, Ahmadi and Arani2008; Picchi & Battiato Reference Picchi and Battiato2018) and serves as a useful archetype for future studies that would address this issue. Indeed, there is no definite model that allows one to take into account the physics in the contact line neighbourhood, and this deserves a specific analysis that is beyond the scope of the present work.

Assuming that no specific viscosity effects occur at the fluid–fluid interface, the stress jump at $\mathscr {A}_{\beta \gamma }$ is compensated by capillary forces, leading to the following boundary condition:

(2.8c)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{p_\beta}-{\boldsymbol{\mathsf{T}}}_{p_\gamma})= \nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}.\end{equation}

After application of Gray's decomposition to the pressure in each phase, the above equation can be equivalently written as

(2.8d)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{\tilde{p}_\beta}-{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}) = \boldsymbol{n}_{\beta\gamma} (\langle p_\beta\rangle^\beta -\langle p_\gamma\rangle^\gamma )_{\boldsymbol{r}_{\beta\gamma}}+ \nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}.\end{equation}

In the above equations, $\boldsymbol {n}_{\beta \gamma }$ is the unit normal vector at $\mathscr {A}_{\beta \gamma }$, directed from the $\beta$ phase towards the $\gamma$ phase, $\boldsymbol {r}_{\beta \gamma }$ locates points at $\mathscr {A}_{\beta \gamma }$ with respect to a fixed system of coordinates (see figure 1b), $H$ is the mean curvature of the interface and $\gamma$ is the interfacial tension between the two fluids. In this boundary condition, possible spatial variations of $\gamma$ along $\mathscr {A}_{\beta \gamma }$ are considered, expressed by $\nabla _s \gamma$, where $\nabla _s$ is the surface gradient operator defined as $\nabla _s=({\boldsymbol{\mathsf{I}}}-\boldsymbol {n}_{\beta \gamma } \boldsymbol {n}_{\beta \gamma })\boldsymbol{\cdot} \boldsymbol{\nabla}$. Note that $2H=-\nabla _s\boldsymbol{\cdot} \boldsymbol {n}_{\beta \gamma }$.

Using Taylor series expansions for the average pressure terms (see (2.4)) about the phase barycentres, $\boldsymbol {x}_\alpha$, recalling that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ can be assumed constant within the unit cell, the boundary condition given in (2.8d) can be rewritten as

(2.8e)\begin{align} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{\tilde{p}_\beta}-{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}) &= \boldsymbol{n}_{\beta\gamma} ( \langle p_\beta\rangle^\beta|_{\boldsymbol{x}_\beta} -\langle p_\gamma\rangle^\gamma|_{\boldsymbol{x}_\gamma} )-\boldsymbol{n}_{\beta\gamma} (\boldsymbol{x}_{\beta}-\boldsymbol{x}_{\gamma})\boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\beta \rangle^\beta\nonumber\\ &\quad +\nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}. \end{align}

While writing (2.8e), use was made of the fact that, if a periodic unit cell representative of the local process is employed for the averaging domain, the macroscopic pressure gradients are equal (see the proof in Appendix B in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022)).

The problem statement in the unit cell is completed with the following periodicity conditions:

(2.8f)\begin{equation} \psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = \boldsymbol{v}_{\alpha}, \tilde{p}_{\alpha}, \end{equation}

and with the average constraint

(2.8g)\begin{equation} \langle \tilde{p}_\alpha\rangle^\alpha=0,\end{equation}

which results from applying the intrinsic averaging operator to (2.3) (with $\psi _\alpha = p_\alpha$) and the use of the Taylor series expansion (2.4), while recalling that $\langle \boldsymbol {z}_\alpha \rangle ^\alpha =\boldsymbol {0}$ and that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\beta$ is taken as a constant in the unit cell.

To finalise this section, it is worth mentioning that the assumptions imposed in this work do not constitute a severe restriction in terms of application of the model. In fact, some systems where it is reasonable to adopt these assumptions (and hence where the macroscopic models presented below are applicable) are the following: incompressible gas–liquid and liquid–liquid packed columns for extraction processes and distillation (Chhabra & Basavaraj Reference Chhabra and Basavaraj2019), incompressible gas and liquid flow through porous shale reservoirs (Singh & Cai Reference Singh and Cai2019), water and non-Newtonian soil flooding for enhanced oil recovery (Sorbie Reference Sorbie1991; Nilsson et al. Reference Nilsson, Kulkarni, Gerberich, Hammond, Singh, Baumhoff and Rothstein2013), as well as industrial applications in packed bed reactors (Giri & Majumder Reference Giri and Majumder2014) and bioreactors (Sen et al. Reference Sen, Nath and Bhattacharjee2017).

3.1. Averaged mass equations

Derivation of the macroscopic mass equation in each phase was reported in previous works dedicated to two-phase flow in homogeneous porous media (Whitaker Reference Whitaker1986, Reference Whitaker1994; Lasseux et al. Reference Lasseux, Quintard and Whitaker1996; Lasseux & Valdés-Parada Reference Lasseux and Valdés-Parada2022) and remains unchanged in the present configuration. The upscaled equation can be written as

(3.1)\begin{equation} \frac{\partial\varepsilon_\alpha}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{v}_\alpha\rangle_\alpha=0.\end{equation}

This expression, which is obtained without assuming that $\mathscr {V}$ is periodic, takes into account the temporal changes of the $\alpha$ phase volume fraction. Time dependency at the macroscale arises from the possible non-stationary character of $\mathscr {A}_{\beta \gamma }$, implicitly making the pore-scale problem time-dependent, even if temporal acceleration is neglected with respect to viscous effects, as already discussed in § 2. Indeed, if $\mathscr {V}$ is assumed periodic, $({1}/{V})\int _{\mathscr {A}_{\beta \gamma }}\boldsymbol {n}_{\beta \gamma } \boldsymbol{\cdot} \boldsymbol {v}_\alpha \, {\rm d} A=\langle \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol {v}_\alpha \rangle _\alpha$, which, after application of the superficial averaging operator to (2.1a), allows concluding that $\boldsymbol{\nabla} \boldsymbol{\cdot} \langle \boldsymbol {v}_\alpha \rangle _\alpha =0$. An equivalent result is obtained if the fluid–fluid interface is stationary, which means

(3.2)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w}=0, \quad \text{at } \mathscr{A}_{\beta\gamma}. \end{equation}

In that case, it follows from the general transport theorem that ${\partial \varepsilon _\alpha }/{\partial t}=0$. This allows concluding that

(3.3)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{v}_\alpha\rangle_\alpha=0, \quad \mathscr{V}\ \text{periodic or } \mathscr{A}_{\beta\gamma} \ \text{stationary}. \end{equation}

3.2. Averaged momentum equations

To obtain the macroscopic momentum equation in each phase, the upscaling approach detailed in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) is used. In brief, it consists of the following steps. First, propose pertinent adjoint (or closure) problems in a periodic unit cell, which can be subsequently related to the flow problem by means of Green's formula. Second, substitute the corresponding differential equations and boundary conditions to carry out the necessary algebraic manipulations to obtain the upscaled models. In this way, the derivations commence by introducing the four second-order tensors, ${\boldsymbol{\mathsf{D}}}_{\alpha \kappa }$, and vectors, $\boldsymbol {d}_{\alpha \kappa }(\alpha,\kappa=\beta,\gamma)$, that solve the following adjoint closure problems in a periodic unit cell (cf. Bottaro Reference Bottaro2019):

Problem I

(3.4a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha \beta}=\boldsymbol{0}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.4b)$$\begin{gather}-\frac{\rho_\alpha}{\mu_{{ref}\alpha}}\boldsymbol{v}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha \beta}=\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\alpha\beta}} +\delta^K_{\alpha\beta} {\boldsymbol{\mathsf{I}}}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.4c)$$\begin{gather}\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} \mu_{{ref}\beta} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\beta\beta}} =\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\gamma} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\gamma\beta}} , \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.4d)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{\beta\beta}={\boldsymbol{\mathsf{D}}}_{\gamma\beta}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.4e)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{ \alpha\beta}=\boldsymbol{0}, \quad\text{at }\mathscr{A}_{ \alpha\sigma}, \end{gather}$$

(3.4f)$$\begin{gather}\psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = {\boldsymbol{\mathsf{D}}}_{\alpha\beta}, \boldsymbol{d}_{\alpha\beta}, \end{gather}$$

(3.4g)$$\begin{gather}\boldsymbol{d}_{\alpha\beta} =\boldsymbol{0}, \quad \text{at }\boldsymbol{r}_{\alpha}^0. \end{gather}$$

Problem II

(3.5a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\gamma}=\boldsymbol{0}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.5b)$$\begin{gather}-\frac{\rho_\alpha}{\mu_{{ref}\alpha}}\boldsymbol{v}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha\gamma}=\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\alpha\gamma}}+ \delta^K_{\alpha\gamma}{\boldsymbol{\mathsf{I}}} , \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.5c)$$\begin{gather}\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\beta} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\beta\gamma}} =\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\gamma} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\gamma\gamma}}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.5d)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{\beta\gamma}={\boldsymbol{\mathsf{D}}}_{\gamma\gamma}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.5e)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{ \alpha\gamma}=\boldsymbol{0}, \quad\text{at }\mathscr{A}_{ \alpha\sigma}, \end{gather}$$

(3.5f)$$\begin{gather}\psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = {\boldsymbol{\mathsf{D}}}_{\alpha\gamma},\boldsymbol{d}_{\alpha\gamma}, \end{gather}$$

(3.5g)$$\begin{gather}\boldsymbol{d}_{\alpha\gamma} =\boldsymbol{0}, \quad \text{at }\boldsymbol{r}_{\alpha}^0. \end{gather}$$

In the above equations, $\boldsymbol {r}_{\alpha }^0$ locates an arbitrary point in the $\alpha$ phase and $\delta _{\alpha \kappa }^K$ is the Kronecker delta. In addition,

(3.6)\begin{equation} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\alpha\kappa}}={-}{\boldsymbol{\mathsf{I}}}\boldsymbol{d}_{\alpha\kappa}+ \frac{\mu(\varGamma_\alpha)}{\mu_{{ref}\alpha}} (\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha\kappa}+\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha \kappa}^{T1}), \end{equation}

where the superscript $T1$ denotes the transpose of a third-order tensor that permutes its first two indices, i.e. $(\boldsymbol{\nabla} {\boldsymbol{\mathsf{D}}}_{\alpha \kappa })_{ijk}^{T1} =(\boldsymbol{\nabla} {\boldsymbol{\mathsf{D}}}_{\alpha \kappa })_{jik}$. Note that the boundary conditions given in (3.4g) and (3.5g) are necessary in order to well pose the closure problems. A practical solution of closure problems I and II, written under the forms given in (3.4) and (3.5), entails a pre-determination of the pore-scale velocity and viscosity fields over the corresponding periodic domain. These adjoint problems can be derived from the associated velocity Green's function pair problems, as detailed in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) for Newtonian creeping two-phase flow. It is worth noting that, in that case, the same closure problems as those previously reported in Auriault & Sanchez-Palencia (Reference Auriault and Sanchez-Palencia1986), Whitaker (Reference Whitaker1994), Lasseux et al. (Reference Lasseux, Quintard and Whitaker1996) and Picchi & Battiato (Reference Picchi and Battiato2018) were retrieved, albeit the macroscale model is different. Analogously, the closure problems reported in (3.4) and (3.5) are the adjoint versions of those derived in Lasseux et al. (Reference Lasseux, Ahmadi and Arani2008) (see (114) and (115) therein) in the case of Newtonian and inertial two-phase flow in homogeneous porous media.

The next step in the derivations consists of relating the associated closure problems with the flow problem with the help of a Green's formula. Considering a constant scalar $d$ and arbitrary and sufficiently regular scalar fields $a$ and $c$, vector fields $\boldsymbol {a}$, $\boldsymbol {b}$ and $\boldsymbol {c}$ and a second-order tensor field ${\boldsymbol{\mathsf{B}}}$, together with the following definitions:

(3.7a)$$\begin{gather} {\boldsymbol{\mathsf{T}}}_{a} ={-}{\boldsymbol{\mathsf{I}}}a + c (\boldsymbol{\nabla} \boldsymbol{a}+ \boldsymbol{\nabla} \boldsymbol{a}^T), \end{gather}$$

(3.7b)$$\begin{gather}{\boldsymbol{\mathsf{T}}}_{{\boldsymbol{b}}} ={-}{\boldsymbol{\mathsf{I}}}\boldsymbol{b}+ c(\boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}}+\boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}}^{T1}), \end{gather}$$

this formula is written (see the proof in Appendix A in Sánchez-Vargas, Valdés-Parada & Lasseux (Reference Sánchez-Vargas, Valdés-Parada and Lasseux2022))

(3.8)\begin{align} &\int_{\mathscr{V}_\alpha } [ \boldsymbol{a} \boldsymbol{\cdot} ( d \boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}} +\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{b}} ) - ({-}d \boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{a} + \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{a}) \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} -\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{a} \boldsymbol{b} \nonumber\\ &\quad +d (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{c} \boldsymbol{a})\boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} +a\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} ]\, {\rm d} V\nonumber\\ &\qquad = \int_{\mathscr{A}_\alpha } [\boldsymbol{a} \boldsymbol{\cdot} ( \boldsymbol{n}_\alpha\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{b}} )- \boldsymbol{n}_\alpha \boldsymbol{\cdot} ({-}d \boldsymbol{c}\boldsymbol{a} +{\boldsymbol{\mathsf{T}}}_{{a}} ) \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}}]\, {\rm d} A. \end{align}

Here, $\mathscr {A}_\alpha$ denotes the surfaces bounding $\mathscr {V}_\alpha$.

Attention is now focused on the derivation of the average momentum balance in the $\beta$ phase. To this end, the above Green's formula is employed, taking $\boldsymbol {a}=\boldsymbol {v}_\beta$, $a= \tilde {p}_\beta$, ${\boldsymbol{\mathsf{B}}}={\boldsymbol{\mathsf{D}}}_{\beta \beta }$, $\boldsymbol {b}=\mu _{{ref}\beta }\boldsymbol {d}_{\beta \beta }$, $c=\mu (\varGamma _\beta )$, $\boldsymbol {c}=\boldsymbol {v}_\beta$ and $d=\rho _\beta$. Once the corresponding differential equations and boundary conditions from the pore-scale problem and closure problem I are substituted in the ensuing relationship, keeping in mind the definition of the superficial averaging operator and the fact that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ and $\rho _\alpha \boldsymbol {b}_\alpha$ are constant in $\mathscr {V}_\alpha$, the following expression for $\langle \boldsymbol {v}_\beta \rangle _\beta$ is obtained:

(3.9a)\begin{align} \langle\boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_{\beta} + \frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta}\, {\rm d} A \nonumber\\ &\quad -\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} )-\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\beta\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma}}\boldsymbol{n}_{\beta \gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta }\, {\rm d} A (\langle p_\beta \rangle^\beta |_{\boldsymbol{x}_\beta} -\langle p_\gamma \rangle^\gamma |_{\boldsymbol{x}_\gamma} -(\boldsymbol{x}_\beta -\boldsymbol{x}_\gamma) \boldsymbol{\cdot} \boldsymbol{\nabla} \langle p_\beta \rangle^\beta). \end{align}

Since ${\boldsymbol{\mathsf{D}}}_{\beta \beta }$ is a solenoidal field that is periodic and zero at $\mathscr {A}_{\beta \sigma }$, one can make use of the divergence theorem to conclude that $\int _{\mathscr {V}_\beta }\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta \beta }\, {\rm d} V=\int _{\mathscr {A}_{\beta \gamma }}\boldsymbol {n}_{\beta \gamma } \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta \beta }\, {\rm d} A=\boldsymbol {0}$. This leads to simplifying the expression of $\langle \boldsymbol {v}_\beta \rangle _\beta$ into the following form:

(3.9b)\begin{align} \langle\boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_{\beta} +\frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta}\, {\rm d} A \nonumber\\ &\quad -\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} ) -\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\beta\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A. \end{align}

This last relationship is still non-local due to the presence of $\tilde {p}_\gamma$ in the last interfacial integral term. To progress towards a local form, Green's formula (3.8) can be employed once more taking $\boldsymbol {a}=\boldsymbol {v}_\gamma$, $a= \tilde {p}_\gamma$, ${\boldsymbol{\mathsf{B}}}={\boldsymbol{\mathsf{D}}}_{\gamma \beta }$, $\boldsymbol {b}=\mu _{{ref}\gamma }\boldsymbol {d}_{\gamma \beta }$, $c=\mu (\varGamma _\gamma )$, $\boldsymbol {c}=\boldsymbol {v}_\gamma$ and $d=\rho _\gamma$. Using the definition of the superficial averaging operator, the resulting expression can be written as follows:

(3.10)\begin{align} &\frac{1}{V}\int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot}( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} )- \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\gamma\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A \nonumber\\ &\quad =(\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma-\rho_\gamma\boldsymbol{b}_\gamma )\boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\gamma\beta}\rangle_\gamma. \end{align}

Substituting this expression in (3.9b), while recalling the interfacial boundary condition given in (3.4d), yields

(3.11)\begin{align} \langle \boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_\beta -\frac{1}{\mu_{{ref}\beta}}(\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma-\rho_\gamma\boldsymbol{b}_\gamma )\boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\gamma\beta} \rangle_\gamma\nonumber\\ &\quad +\frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta} \, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w} \boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\beta\beta}\,{\rm d} A. \end{align}

The last area integral term on the right-hand side of the above equation is related to the displacement of $\mathscr {A}_{\beta \gamma }$ when inertia is present.

The above procedure can be repeated in a similar way for the $\gamma$ phase. The resulting equation can be written as follows:

(3.12)\begin{align} \langle\boldsymbol{v}_\gamma\rangle_\gamma &={-}\frac{1}{\mu_{{ref}\gamma}}(\boldsymbol{\nabla}\langle p_\gamma\rangle^\gamma -\rho_\gamma\boldsymbol{b}_\gamma ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\rangle_\gamma - \frac{1}{\mu_{{ref}\gamma}}(\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\beta\gamma}\rangle_\beta\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\gamma}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\gamma}V} \int_{\mathscr{A}_{\beta\gamma} }\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot}\boldsymbol{w}\boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\, {\rm d} A. \end{align}

At this point, it is convenient to adopt the following nomenclature ($\alpha,\kappa =\beta, \gamma$):

(3.13)\begin{equation} {\boldsymbol{\mathsf{H}}}_{\alpha\kappa}=\frac{\mu_{{ref} \kappa}}{\mu_{{ref}\alpha}}\langle{\boldsymbol{\mathsf{D}}}_{\kappa \alpha}\rangle_\kappa^T,\end{equation}

${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ and ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$ being respectively identified as the classical dominant and coupling apparent permeability tensors in the $\alpha$ phase. This allows writing the macroscopic momentum balance equations in the following compact form:

(3.14a)\begin{align} \langle\boldsymbol{v}_\alpha\rangle_\alpha &={-} \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\alpha}}{\mu_{{ref}\alpha}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha -\rho_\alpha \boldsymbol{b}_\alpha ) -\frac{{\boldsymbol{\mathsf{H}}}_{\alpha\kappa}}{\mu_{{ref}\kappa}}\boldsymbol{\cdot} (\boldsymbol{\nabla}\langle p_\kappa\rangle^\kappa -\rho_\kappa\boldsymbol{b}_\kappa )\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\alpha}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\alpha}V} \int_{\mathscr{A}_{\beta\gamma} } \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w} \boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A. \end{align}

This is the first important result of the present work. The first two terms on the right-hand side of (3.14a) are the Darcy-like terms. The first one accounts for the viscous resistance in the phase under concern, while the second one results from the viscous coupling between the two fluid phases. These terms include bulk inertial effects, justifying the apparent and non-intrinsic character of the permeability tensors ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$. In addition to the Darcy-like terms, a contribution of the interfacial effects is also present, which is accounted for by the last two terms in (3.14a). The first one, dealing with interfacial tension effects, is formally the same as that identified in a previous work dedicated to creeping Newtonian two-phase flow in porous media (Lasseux & Valdés-Parada Reference Lasseux and Valdés-Parada2022). The second one results from inertial effects. Note that, both inertial and non-Newtonian effects are reflected in the values of the closure variable ${\boldsymbol{\mathsf{D}}}_{\alpha \alpha }$ at $\mathscr {A}_{\beta \gamma }$ (and also in $H$). The first interfacial term was not present in the development of the average model proposed by Lasseux et al. (Reference Lasseux, Ahmadi and Arani2008) in the case of two-phase inertial Newtonian flow. This is because a simplification of the contribution of the capillary forces was assumed. Such an assumption lacks careful justification, and this was extensively discussed in the case of creeping Newtonian two-phase flow (see Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022), pp. 14–15). When the fluid–fluid interface is stationary, or flow remains in the creeping regime in the $\alpha$ phase, this equation reduces to

(3.14b)\begin{align} \langle\boldsymbol{v}_\alpha\rangle_\alpha &={-}\frac{{\boldsymbol{\mathsf{H}}}_{\alpha\alpha}}{\mu_{{ref}\alpha}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha -\rho_\alpha\boldsymbol{b}_\alpha ) - \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\kappa}}{\mu_{{ref}\kappa}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\kappa\rangle^\kappa -\rho_\kappa\boldsymbol{b}_\kappa )\nonumber\\ &\quad +\frac{1}{\mu_{{ref}\alpha}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A,\nonumber\\ &\quad\quad \quad \quad \quad \text{creeping flow in the }\alpha \text{ phase or}\ \mathscr{A}_{\beta\gamma} \ \text{stationary}. \end{align}

Note that (3.14b) is valid in the creeping flow regime, even if $\boldsymbol {n}_{\beta \gamma }\boldsymbol{\cdot} \boldsymbol {w}\ne 0$. To that extent, the model reported in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) is valid even if the fluid–fluid interface is not in its stationary configuration.

All the terms in the macroscopic momentum equations are determined from the solution of the two (adjoint) closure problems I and II, which requires knowledge of the interface position, $\mathscr {A}_{\beta \gamma }$, and of its speed of displacement, $\boldsymbol {w}$, if inertia is significant and the interface has not reached a stationary configuration. Note that these closure problems also require accounting for the velocity field in each phase even under creeping flow conditions for the computation of the viscosity coefficient $\mu (\varGamma _\alpha )$.

3.3. Symmetry and reciprocity relationships for ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ and ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$ ($\alpha, \gamma =\beta, \gamma$, $\alpha \ne \kappa$)

In Appendix A, a symmetry analysis of ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ is provided, along with a reciprocity relationship between ${\boldsymbol{\mathsf{H}}}_{\beta \gamma }$ and ${\boldsymbol{\mathsf{H}}}_{\gamma \beta }$. The results are respectively given in (A7) and (A10). The relationship (A7) shows, first, that ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ is not symmetric and, second, that its skew-symmetric part only results from inertial effects. In addition, the reciprocity relationship given in (A10) indicates that the two coupling apparent permeability tensors are not independent and that the right-hand side also only results from inertia. Indeed, in the absence of inertial effects, the right-hand sides of both (A7) and (A10) are zero. This allows concluding that, under creeping flow conditions, the dominant permeability tensors are symmetric, as for Newtonian flow, and the coupling permeability tensors satisfy the same relationship as for $\mu (\varGamma _\alpha )={Cst}=\mu _{{ref}\alpha }$. This extends the results reported in Lasseux et al. (Reference Lasseux, Quintard and Whitaker1996) and Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017) (section II.E) for Newtonian two-phase flow, showing that the non-Newtonian character of the flow does not alter the properties of the permeability tensors in the absence of inertia. It also extends the result obtained for one-phase incompressible Newtonian and inertial flow in porous media for which the apparent permeability tensor was shown to be non-symmetric, with its skew-symmetric part only due to inertia (Lasseux & Valdés-Parada Reference Lasseux and Valdés-Parada2017, § II.A).

To complete the macroscopic description of the two-phase flow process, it is of interest to derive an expression for the macroscpic pressure difference. This is the objective of the development that follows.

6.1. Analysis of the macroscopic velocity

This part of the work is dedicated first to the study of the influence of rheological parameters under creeping flow conditions for the configuration in which the $\beta$ phase is continuous (see figure 2a). The use of creeping flow conditions in the first set of simulations is justified both to contrast with previously reported results (Lasseux & Valdés-Parada Reference Lasseux and Valdés-Parada2022), and also to isolate the effects of rheology. Secondly, inertial effects are studied in the three configurations depicted in figure 2. In all cases, emphasis is laid upon validation with DNS results.

6.1.1. Influence of the rheological parameters in continuous fluid phases

In figure 3, a first set of results corresponding to the pore-scale flow streamlines, superimposed to the velocity magnitude, are represented, taking $n=0.4$, $\lambda ^*=0.5$ and four values of $S_\beta$ ranging from 0.4 to 0.7. For the Newtonian case (figure 3a,c,e,g), the results reproduce those reported in figure 5 of Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) and are recalled here for comparison purposes. As expected, $\mathscr {A}_{\beta \gamma }$ coincides with a streamline in that case. This feature is not observed in figure 3(b,d, f,h), which correspond to a situation where the $\beta$ phase obeys a Carreau model. In addition, the viscosity parameters employed here, corresponding to a shear-thinning wetting phase, lead to a very significant increase of the velocity by more than an order of magnitude in both phases in comparison with the case in which the $\beta$ phase is Newtonian. The change in the flow structure observed when the $\beta$ phase is non-Newtonian results from both the change in rheology and the fact that the fluid–fluid interface does not correspond to its steady configuration. In all cases, the largest wetting-phase velocity is located in the constriction between the solid cylinder and $\mathscr {A}_{\beta \gamma }$, which corresponds to a low-pressure zone. In contrast, the portions of the wetting phase where eddies are present correspond to those where the viscosity is the largest and correlate to large-pressure zones. Consequently, as the saturation increases, the low-pressure zones increase in size, thus favouring the overall fluid flow. The same effect is observed when $\lambda ^*$ increases, or when $n$ decreases, which is consistent with the shear-thinning nature of the Carreau fluid considered here.

Figure 3. Fields of the velocity magnitude and streamlines (in white) under Newtonian (a,c,e,g) and non-Newtonian (b,d, f,h) flow conditions for the wetting $\beta$ phase, taking four values of the $\beta$-phase saturation ($S_\beta$): (a,b) $S_\beta =0.4$; (c,d) $S_\beta =0.5$; (e, f) $S_\beta =0.6$; (g,h) $S_\beta =0.7$. The fluid–fluid interface is represented as a black curve and the central solid cylinder in white. Simulations for non-Newtonian flow (b,d, f,h) correspond to a Carreau fluid considering the power-law index $n=0.4$ and the dimensionless relaxation time $\lambda ^*=5$. Due to symmetry, results are only reported for the top half of the unit cell depicted in figure 2(a). Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

At this point, it is worth pondering the contribution of the different terms present in the dimensionless macroscopic flow model as expressed in (B7). In table 1, data are reported to evaluate both the Darcy-like part and the interfacial part of $\langle v_{\alpha x}^*\rangle _\alpha$ predicted by this model, taking $n=0.4$ and $\lambda ^*=5$, for the same wetting-phase saturation values considered in figure 3. These predictions yield relative errors with respect to DNS of less than 1 % for all the values of $S_\beta$ considered here. Moreover, the Darcy-like terms are found to be, at least, one order of magnitude larger than the interfacial terms for the set of parameters considered here. Furthermore, since, for each given saturation, the macroscopic pressure gradient is along $\boldsymbol {e}_x$ (see (6.2)), $\langle v_{\alpha y}^*\rangle _\alpha$ is numerically verified to be zero, so, from this point forth, focus is laid upon $\langle v_{\alpha x}^*\rangle _\alpha$. Results in table 1 also show that the contributions from $K^*_{\beta \beta xx}$, $K^*_{\beta \gamma xx}$ and $K^*_{\gamma \beta xx}$ increase with $S_\beta$. However, for $K^*_{\gamma \gamma xx}$ the same observation is not applicable and this needs to be further analysed. To address this point, figures 4 and 5 report the predictions of the permeability coefficients and the average velocities for different wetting-phase saturations, while varying $n$ and $\lambda ^*$. In these figures, results are presented in terms of the $xx$ component $k_{r\alpha \kappa }$ of the relative permeability tensor, defined as

(6.3)\begin{equation} k_{r\alpha\kappa}=\frac{K_{\alpha\kappa xx}}{K}, \quad \alpha,\kappa=\beta,\gamma, \end{equation}

where ${K}_{\alpha \kappa xx}$ is the $xx$ component of the permeability tensor ${\boldsymbol{\mathsf{K}}}_{\alpha \kappa }$ and $K$ is the intrinsic permeability, which is a scalar since, for the case under consideration, the structure is orthotropic.

Table 1. Predictions of the macroscopic flow model for different saturations of the non-Newtonian wetting phase, $S_\beta$. The error is defined as $\text {Error}_\alpha ~(\%)= 100 \times |\langle v_{\alpha x}\rangle _{\alpha DNS}-\langle v_{\alpha x}\rangle _\alpha |/\langle v_{\alpha x}\rangle _{\alpha DNS}$. Here, $K^*_{\alpha \kappa xx}$ stands for the xx component of the corresponding dimensionless permeability tensor and $\alpha _\alpha ^*={2\mu _{{ref}\gamma }}/({\mu _{{ref}\alpha } Ca V^*})\int _{\mathscr {A}_{\beta \gamma } } H^* \boldsymbol {n}_{\beta \gamma }\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}^*_{\alpha \alpha }\, {\rm d} A^*\boldsymbol{\cdot} \boldsymbol {e}_x$. Also, $\varepsilon =0.8$, $Ca=10$, $\mu ^*=0.1$, $n=0.4$ and $\lambda ^*=5$. The results correspond to the geometry depicted in figure 2(a).

Figure 4. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($k_{r\alpha \alpha }$) and coupling ($k_{r\alpha \kappa }$, $\alpha \neq \kappa$) relative permeabilities versus $S_\beta$, taking the intrinsic permeability as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $\lambda ^*=5$. $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$. The results correspond to the geometry depicted in figure 2(a).

Figure 5. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($k_{r\alpha \alpha }$) and coupling ($k_{r\alpha \kappa }$, $\alpha \neq \kappa$) relative permeabilities versus $S_\beta$, taking the intrinsic permeability as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$. $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$; $\bigstar$, $\lambda ^*=100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$. The results correspond to the geometry depicted in figure 2(a).

In figure 4, an analysis is presented fixing $\lambda ^*$ and varying the power-law index $n$ in a range of values representative of different non-Newtonian fluids. Regarding these results, the following remarks are in order.

• • From figure 4(a,b), it is clear that the values of the x component of the velocity, resulting from averaging the pore-scale simulation fields and predicted from solving the upscaled model, show excellent agreement for both fluid phases. In fact, the relative error of the upscaled model predictions with respect to those from pore-scale simulations is less than 5 % in all cases. This comment is also applicable to the results reported in figure 5(a,b), for $0\leq \lambda ^*\leq 100$.

• • Clearly, as the non-Newtonian character of the wetting phase is more pronounced, the average flow in both phases is favoured, in agreement with the shear-thinning model under consideration. The increment in both phases is very important as it can reach two orders of magnitude with respect to the Newtonian case. This is attributed to the increase of the relative permeability xx component and this is in agreement with the observations made in figure 3.

• • As $n$ decreases, the relative permeabilities and the average velocities become more sensitive to its variations. This is more evident for $n<0.5$.

• • The graphs in figure 4(e, f), describing the xx component of the coupling relative permeabilities, are related by $k_{r\gamma \beta }= \mu ^* k_{r\beta \gamma }$. This proportionality can be deduced from (A10) under creeping flow conditions.

The second degree of freedom in this analysis is the dimensionless relaxation time, $\lambda ^*$. Results corresponding to the sensitivity of macroscale fluid velocities and relative permeabilities to this parameter are reported in figure 5 for several values of the wetting-phase saturation. Regarding these results, the following observations are pertinent.

• • Consistent with the results shown in figure 4, flow within the unit cell is improved as the shear-thinning character of the wetting phase is increased. Once more, both the relative permeabilities and macroscale velocities increase by roughly two orders of magnitude as $\lambda ^*$ increases from 1 to 100. Note that the flow is more sensitive to variations of $\lambda ^*$ in the range $0\le \lambda ^*\le 5$.

• • Similar to the analysis of the impact of $n$, it must be noted that $k_{r\beta \gamma }/k_{r\gamma \beta }=10$ (see figure 5e, f) and the similarity with figure 4(e, f). This confirms the expected reciprocity relationship reported in § 3.3 (see also (A10)), in agreement with the fact that $\mu ^*=0.1$.

The above analysis shows the important role played by both the dominant and coupling Darcy-like terms. For the specific combination of values of the capillary number and viscosity ratio, together with the porous medium geometry and interface configuration chosen here, there is no relevant contribution from surface tension effects (i.e. $\boldsymbol {\alpha }^*_\alpha$, cf. (B8a)) to the prediction of the macroscopic fluid velocities. However, this should not be considered as a general result as the contribution of this interfacial term was shown to become even dominant in some circumstances, in particular while decreasing $Ca$ or increasing $\mu ^*$ (see table 1 in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022)). The contribution of this term, along with the inertial interfacial term, $\boldsymbol {\omega }^*_\alpha$ (cf. (B8)), is further explored next when inertia is present.

6.1.2. Influence of inertial effects

Let now attention be directed to the influence of inertial effects over the predictions of the macroscale velocity and effective-medium coefficients. The results presented below require specifying the Reynolds number values, $Re_\alpha$, $\alpha =\beta,\gamma$, in both fluid phases. For this reason, $Re_\beta$ was arbitrarily set equal to one-tenth of $Re_\gamma$. Despite the fact that in the case of figure 2(c) the absence of symmetry of the unit cell makes the apparent permeability tensors fully populated, in the following, only the $xx$ components of these tensors are reported, normalised with the $xx$ component of the intrinsic permeability. This is,

(6.4)\begin{equation} h_{r\alpha \kappa} = \frac{H_{\alpha \kappa xx}}{K}, \quad \alpha, \kappa = \beta, \gamma. \end{equation}

Accordingly, only the $x$ components of the average velocities are reported.

In figures 6 and 7, the predictions of the average velocities and relative permeabilities versus $S_\beta$, taking four values of $Re_\gamma$, resulting from solving the closure problems in the geometries depicted in figures 2(a) and 2(b), are respectively reported for $n=0.5$ and $\lambda ^*=5$ in the Carreau model. From these results, it is clear that the influence of the rheology, combined with inertia, yields a complex dependence of the average velocities and effective-medium coefficients on $S_\beta$ for the Reynolds number values considered here. An important remark must be made for the situation in which the $\beta$ phase is trapped (see figure 2b). It is worth noticing in that case, that $\langle v^*_{\beta x}\rangle _\beta$ is not zero, despite the fact that it is trapped. However, this is not surprising as this is a signature of the fact that the interface is not at its steady position (i.e. the interface position can be thought of as being captured at a given time at which this phase is rearranging, before equilibrium is reached). Therefore, there is a net local momentum transport in the $\beta$ phase in the $x$ direction. This should not be confused with the fact that, in this configuration, the $\beta$ phase remains macroscopically trapped in the porous medium, hence yielding no mass flow rate at the macroscopic system boundaries.

Figure 6. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $S_\beta$, taking the $xx$ component of the intrinsic permeability tensor as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(a) in which the $\beta$ phase is continuous. $\bullet$ (red), creeping flow regime ($Re_\gamma =0$); $\blacktriangle$, $Re_\gamma =10$; $\blacktriangledown$, $Re_\gamma =50$; $\blacklozenge$ (grey), $Re_\gamma =100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

Figure 7. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $S_\beta$, taking the $xx$ component of the intrinsic permeability tensor as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(b) in which the $\beta$ phase is discontinuous. $\bullet$ (red), creeping flow regime ($Re_\gamma =0$); $\blacktriangle$, $Re_\gamma =10$; $\blacktriangledown$, $Re_\gamma =50$; $\blacklozenge$ (grey), $Re_\gamma =100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

The effects of $\alpha ^*_{\kappa x}$ are not reported here since they are negligible when the $\beta$ phase is continuous (cf. figure 2a), and they are exactly zero when it is discontinuous (cf. figure 2b), since $H$ is constant in that case. Note that the dependence of the average velocities and $h^*_{r\alpha \alpha }$ ($\alpha =\beta,\gamma$) on the Reynolds number is not monotonic over the reported range of saturations. This can be attributed to pore-scale balance between viscous dissipation and inertial acceleration in both fluid phases and to the fact that the fluid–fluid interface is not at its stationary position. Nevertheless, as expected, when the $\beta$ phase is continuous, the values of both the average velocities and $h^*_{r\alpha \alpha }$ are larger than those corresponding to the case in which the $\beta$ phase is trapped. This contrast is more pronounced in the $\gamma$ phase, as a result of the fact that, when the $\beta$ phase is trapped, the wetting phase exerts an extra drag on the surrounding non-wetting phase flow. Finally, for $h^*_{r\alpha \kappa }$ ($\alpha, \kappa =\beta,\gamma$, $\alpha \ne \kappa$), the simple reciprocity relationship observed in the previous section is no longer valid, especially when the $\beta$ phase is discontinuous (see figure 7e, f). This is because the terms on the right-hand side of (A10) are no longer vanishingly small compared with their counterparts on the left-hand side when inertia is present. Despite this, the predictions for $h_{r\beta \gamma }$ remain one order of magnitude larger than those for $h_{r\gamma \beta }$ as was observed in the previous paragraphs under creeping flow conditions.

The last part of the analysis of inertial effects is focused on the configuration depicted in figure 2(c). In this case, $\varepsilon =0.915$, $S_\beta =0.346$, $\mu ^*=0.1$, $n=0.5$, $\lambda ^*=5$ and $Ca=10$, so the results for the velocity and the relative (dominant and coupling) apparent permeabilities reported in figure 8 are presented as functions of the Reynolds number in the $\gamma$ phase. These results show that larger values of $Re_\gamma$ are required to appreciate the effect of inertia, in contrast to the results for the configurations in figures 2(a) and 2(b). This is attributed to the viscous drag exerted by the solid on both fluid phases, and this was not the case in the other geometries for which the $\gamma$ phase experienced lubrication effects. The predictions for $\langle v_{\alpha x}^*\rangle _\alpha$ are qualitatively alike in both fluid phases, although $\langle v_{\gamma x}^*\rangle _\gamma$ and $h_{r\gamma \gamma }$ are one order of magnitude larger than their counterparts in the $\beta$ phase. The dependence of $\langle v_{\alpha x}^*\rangle _\alpha$ on the Reynolds number is consistent with observations made in inertial one-phase flow in porous media (cf. Lasseux, Arani & Ahmadi Reference Lasseux, Arani and Ahmadi2011). Finally, the coupling permeability terms maintain the same properties observed in the previous geometries, i.e. $h_{r\beta \gamma }$ is one order of magnitude larger than $h_{r\gamma \beta }$, albeit their dependence upon the Reynolds number is not the same. It is worth adding that the results presented in figures 6–8 show an excellent agreement between DNS and the upscaled model predictions for the average fluid velocities, thus validating the upscaled model in the presence of inertia.

Figure 8. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus $Re_\gamma$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $Re_\gamma$, taking the $xx$ component of the intrinsic permeability as the reference. Results correspond to $Ca=10$ and $\mu ^*=0.1$ for a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(c) in which the $\beta$ phase is both continuous and discontinuous and both phases are in contact with the solid phase $\sigma$. In (a,b), results from DNS and the upscaled model are reported with open and filled circles, respectively. Other parameters are $\varepsilon=0.915$ and $S_\beta=0.346$.

To conclude this part of the analysis, it is pertinent to examine the contributions to $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) from the $x$ component of the interfacial terms $\boldsymbol {\alpha }_\kappa ^*$ and $\boldsymbol {\omega }_\kappa ^*$ ($\kappa =\beta,\gamma$). In figure 9, these contributions are reported normalised with $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) in both fluid phases against $Re_\gamma$. These results show that, for the particular configuration under study, the contributions from $\boldsymbol {\omega }_\kappa ^*$ can be neglected for $Re_\gamma \le 10$ with respect to those from $\boldsymbol {\alpha }_\kappa ^*$. Nevertheless, $\boldsymbol {\omega }_\kappa ^*$ becomes especially relevant for $Re_\gamma >10^3$ as it can be roughly equal to the macroscopic velocity values of both fluid phases. In contrast, the contributions from $\boldsymbol {\alpha }_\kappa ^*$, representing between 3 % and 10 % of the average velocities, remain relevant in the entire range of Reynolds number values considered here.

Figure 9. Per cent contribution to $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) of the interfacial term $x$ component defined in (B8) in (a) the $\beta$ phase and (b) the $\gamma$ phase (i.e. $\psi _\kappa ^{\%} = 100 \times | \psi _{\kappa x}^*|/\langle v_{\kappa x}^*\rangle _\kappa$, $\psi =\alpha,\omega$, $\kappa = \beta,\gamma$) as functions of the Reynolds number in the $\gamma$ phase ($Re_\gamma$). Results correspond to $Ca=10$ and $\mu ^*=0.1$ for a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(c). Other parameters are $\varepsilon=0.915$ and $S_\beta=0.346$.

6.2. Analysis of the average pressure difference

An analysis similar to that reported above for the average velocities is now performed for the macroscopic pressure difference model, as given in (B10). For the sake of brevity, only the case depicted in figure 2(a) in the creeping regime is considered, keeping $\nabla ^*\langle p^*_\alpha \rangle ^\alpha$ as given by (B2). In figure 10(a), the predictions of ${\rm \Delta} \mathcal {P}^*$ from the dimensionless upscaled equation (B10) in which Reynolds numbers are zero are compared with the DNS results, ${\rm \Delta} \mathcal {P}^*_{DNS}$, for different wetting-phase saturations $S_\beta$ and values of the power-law index $n$ ranging from 0.4 to 1.0, fixing $\lambda ^*=5$. Clearly, the derived model predictions perfectly match the pore-scale numerical simulations. In all cases, the error %, calculated as $100 \times |{\rm \Delta} \mathcal {P}^*_{DNS}-{\rm \Delta} \mathcal {P}^*|/ {\rm \Delta} \mathcal {P}^*_{DNS}$, remains smaller than 5 %. This is also true for the results reported in figure 10(b) corresponding to values of the dimensionless relaxation time $\lambda ^*$ ranging from 0 to 50, keeping $n=0.5$. Results for ${\rm \Delta} \mathcal {P}^*$ for larger values of $\lambda ^*$ do not significantly change, and, therefore, are not reported in figure 10(b).

Figure 10. Dimensionless macroscopic pressure difference ${\rm \Delta} \mathcal {P}^*$ versus the wetting-phase saturation $S_\beta$ resulting from the upscaled model (filled symbols) and from DNS (open circles). Results correspond to the geometry depicted in figure 2(a) with a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase. In (a), $\lambda ^*=5$ and $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (b) $n=0.5$ and $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1.$

Increasing the shear-thinning effect of the $\beta$ phase, by either decreasing $n$ or increasing $\lambda ^*$, yields a decrease of the macroscopic pressure difference as shown in figures 10(a) and 10(b), respectively. This is particularly noticeable for $S_\beta <0.6$, whereas, for larger values of the saturation, the impact of the wetting-phase rheology is almost insignificant. Moreover, in agreement with the velocity fields reported in figure 3, it is verified that the average pressure in each phase decreases while increasing the wetting-phase saturation. This is consistent with the fact that the size of the recirculation zones in the wetting phase decreases with $S_\beta$. Nevertheless, the above does not necessarily explain why the macroscopic pressure difference tends to zero as the wetting-phase saturation increases. To address this issue, it is worth recalling that ${\rm \Delta} \mathcal {P}^*$ results from the contribution of both dynamic (i.e. the macroscopic forcing) and capillary effects, represented by $\psi ^*$ and $s_{\beta \gamma }^*$, respectively. Therefore, the impact of the rheological parameters $n$ and $\lambda ^*$ on $\psi ^*$ and $s_{\beta \gamma }^*$, representing the (dimensionless) dynamic and capillary contributions, respectively given by (B11a) and (B11b), is now analysed.

The effective quantities $\psi ^*$ and $s_{\beta \gamma }^*$, obtained from the closure problem solution with $n$ ranging from 1 to 0.4, while fixing $\lambda ^*=5$, are reported versus $S_\beta$ in figures 11(a) and 11(b), whereas, in figures 11(c) and 11(d), the same representation is provided for results obtained by varying $\lambda ^*$ from 0 to 50, fixing $n=0.5$. From figures 11(a) and 11(c), it can be noted that the contribution of the dynamic effects $\psi ^*$ to ${\rm \Delta} \mathcal {P}^*$ is always negative, contrary to the dimensionless capillary effect contribution $s_{\beta \gamma }^*$ (see figures 11b and 11d). From these figures, it is clear that both $\psi ^*$ and $s_{\beta \gamma }^*$ tend to zero as $S_\beta$ increases. This results in a macroscopic pressure difference that also tends to zero when $S_\beta$ takes large enough values (i.e. $>0.6$). In addition, results in figures 11(a) and 11(c) show that $\psi ^*$ increases in magnitude when the shear-thinning effect in the $\beta$ phase becomes more pronounced, whereas the opposite holds for $s_{\beta \gamma }^*$ (see figures 11b and 11d). It must be noted that both dynamic and capillary effects play a significant role in the prediction of the average pressure difference for the conditions under consideration here. This is in contrast with what was concluded for the average velocities for which the interfacial term is insignificant under the same conditions.

Figure 11. Dynamic (i.e. macroscopic forcing) $\psi ^*=\varphi ^*_{\beta x}+\varphi ^*_{\gamma x}+x^*_\beta -x^*_\gamma$ and capillary effect $s_{\beta \gamma }^*$ contributions versus the wetting-phase saturation $S_\beta$. In (a,b) $\lambda ^*=5$ and $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (c,d) $n=0.5$ and $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$. In all cases, $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$ in the geometrical configuration of figure 2(a).

To conclude this section, it is worth remarking that the main objective of the present work is the derivation of the macroscopic models summarised in § 5. Therefore, an exhaustive analysis of all the degrees of freedom involved in the closure problems solution was not carried out. Indeed, the numerical analysis presented here is only illustrative, as the upscaled models are applicable under broader flow conditions. Moreover, rather than reproducing the conditions of a particular experimental system, the interest lies here in analysing the contributions of each term constituting the upscaled models. Importantly, in all the situations considered in this work, the upscaled models for both the average velocities and average pressure difference are validated by the pore-scale direct numerical simulations. It must be emphasised that this holds even if the fluid–fluid interface is not in its stationary configuration, if it is continuous or not through the unit cell and in potentially mixed-wet conditions. This should be regarded as a first validation step and motivates further comparisons with laboratory experiments.