1. Introduction
Modelling two-phase flow in porous media has been a long-lasting subject of research, and remains a challenging task. Over the past two decades, developments in image acquisition and processing (see e.g. Singh et al. Reference Singh, Jung, Brinkmann and Seemann2019, and references therein) as well as progress in direct numerical simulations (DNS) at the pore scale (e.g. Konangi et al. Reference Konangi, Palakurthi, Karadimitriou, Comer and Ghia2021) have been proved to be useful in providing valuable pore-scale information that confirms, and sometimes questions, existing theoretical macroscopic models. Nevertheless, mathematical modelling of two-phase flow in porous media remains of major importance.
A traditional approach for slow, incompressible and immiscible two-phase flow, originating from the pioneering work of Muskat & Meres (Reference Muskat and Meres1936), consists in using the so-called empirical generalized Darcy's equations for momentum transport, together with the macroscopic mass conservation equations. However, to close the macroscopic model, an additional relationship is required (Chavent & Jaffré Reference Chavent and Jaffré1986). On the basis of the assumption that the flow is slow enough to be considered as a succession of equilibria, this closure relationship is taken as the macroscopic pressure difference, often called the macroscopic capillary pressure. This terminology was adopted due to the quasi-equilibrium assumption, so that the pore-scale pressure jump at the fluid–fluid interface is compensated only by capillary effects together with the fact that the pressure in each phase is quasi-uniform and is thus assumed to be equal to that at the fluid–fluid interface. In addition, since the macroscopic flow description is given in terms of the local (wetting-phase) saturation, the macroscopic capillary pressure has been considered as a function of this unique variable.
The shortcomings of such a heuristic description have been widely discussed in the literature over the past four decades. First, the fact that, due to dynamic effects, the difference between the volume-averaged pressures should be distinguished from the interfacial average of the equilibrium pressure difference has been pointed out. In particular, Gray, Bruning & Miller (Reference Gray, Bruning and Miller2019) emphasized that, under dynamic conditions, the macroscopic capillary pressure should be evaluated in terms of the interfacial average of the pressure difference between the two fluids, although, at equilibrium, it corresponds to the volume-average pressure difference. This point of view was shared by Starnoni & Pokrajac (Reference Starnoni and Pokrajac2020), who derived a macroscale Laplace-type equation using volume averaging. Their result is similar to that anticipated by Whitaker (Reference Whitaker1994). Second, the unique dependence of the (equilibrium) macroscopic capillary pressure on the saturation has also been questioned.
On a thermodynamic basis, Hassanizadeh & Gray (Reference Hassanizadeh and Gray1990, Reference Hassanizadeh and Gray1993) developed a relationship for the macroscopic pressure difference that differs from the macroscopic capillary pressure by a dynamic term. Assuming a linear theory, this term was proposed to be proportional to the time rate of variation of the saturation as anticipated by Marle (Reference Marle1981). However, a closure relationship to estimate the corresponding proportionality coefficient was not provided. Furthermore, these authors pinpointed the dependence of the (equilibrium) macroscopic capillary pressure on many parameters, which, under classical flow conditions, can be restricted to the fluid–fluid interfacial area and the saturation. This allowed them to conclude that the equilibrium macroscopic capillary pressure should be viewed as a unique surface parametrized by the interfacial specific area and saturation, and that the hysteretic dependence of the capillary pressure on the saturation, which has been widely reported, results from a projection of this surface onto the saturation axis.
Experimental measurements in transparent micromodels of the equilibrium macroscopic pressure were carried out by Pyrak-Nolte et al. (Reference Pyrak-Nolte, Nolte, Chen and Giordano2008) to investigate functionalities proposed in the above two references. Their results highlighted the fact that the terms of the capillary pressure related to the interfacial area were nearly constant for any equilibrium situation, and, consequently, that terms associated with the phase-free energy variations with respect to the saturation and related to saturation gradients play a very important role. Under dynamic conditions, Joekar-Niasar, Hassanizadeh & Dahle (Reference Joekar-Niasar, Hassanizadeh and Dahle2010) and Joekar-Niasar & Hassanizadeh (Reference Joekar-Niasar and Hassanizadeh2012) used a three-dimensional pore-network modelling on a regular lattice and confirmed that a classical capillary pressure function of the saturation was inappropriate. They concluded that a unique surface relating the capillary pressure to the saturation and interfacial area was in fair agreement with the observations, although the dependence upon the capillary number was not fully captured.
Using the thermodynamically constrained averaging theory, Gray & Miller (Reference Gray and Miller2011) proposed an expression for the capillary pressure evolution to equilibrium that incorporates the rate of change of the specific fluid–fluid interfacial area. Since no closure scheme was provided to determine the weighting coefficients of both the rate of change of saturation and interfacial specific area, their determination was left to experimental investigation.
In a recent development proposed by Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022), a formal derivation of the macroscopic momentum and mass conservation equations for isothermal, creeping Newtonian flow of two immiscible and incompressible phases, 
$\beta$ and 
$\gamma$, in homogeneous porous media, was obtained, which reads (see the definitions of the different quantities in § 2)
\begin{gather} \langle {\boldsymbol{v}}_\alpha\rangle_\alpha ={-}\frac{{\boldsymbol{\mathsf{K}}}_{\alpha\alpha}^*}{\mu_\alpha} \boldsymbol{\cdot} \boldsymbol{\nabla} \mathscr{P}_\alpha -\frac{{\boldsymbol{\mathsf{K}}}_{\alpha\kappa}^*}{\mu_\kappa} \boldsymbol{\cdot} \boldsymbol{\nabla} \mathscr{P}_\kappa +\frac{1}{\mu_\alpha V} \int_{\mathscr{A}_{\beta \gamma}} (2\gamma H {\boldsymbol{n}}_{\beta\gamma} + \boldsymbol{\nabla}_s\gamma) \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\,{\rm d}A, \end{gather}
\begin{gather}\frac{\partial \varepsilon_\alpha}{\partial t} +\boldsymbol{\nabla}\boldsymbol{\cdot} \langle\boldsymbol{v}_\alpha\rangle_\alpha=0, \end{gather}
with 
$\alpha,\kappa =\beta,\gamma$ and 
$\kappa \ne \alpha$. In the above momentum equation, the (dominant and coupling) permeability tensors are defined as 
${\boldsymbol{\mathsf{K}}}_{\alpha \kappa }^* = \langle {\boldsymbol{\mathsf{D}}}_{\alpha \kappa }\rangle _\alpha$, where 
${\boldsymbol{\mathsf{D}}}_{\alpha \kappa }$ are closure variables that solve two ancillary closure problems (see (5.14) and (5.15) in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022)).
For completeness, it is of major interest to derive a closed upscaled equation for the average pressure difference between the two fluid phases (which, for the sake of semantic exactness, is not referred to as the dynamic capillary pressure). This is the purpose of the present companion work developed in the next sections as follows. In § 2, the pore-scale flow problem is formulated for a two-phase, Newtonian, incompressible, isothermal and creeping flow. Assuming the existence of an elementary volume representative of the process, the formulation is provided in the corresponding periodic unit cell. Using a simplified version of the volume-averaging method and elements of the adjoint method, an adjoint problem is proposed that is further combined with the physical pore-scale problem with the aid of Green's formula. The procedure, which is free of any simplifying assumption and yields a closed equation for the average pressure difference, is detailed in § 3. Section 4 is dedicated to first evaluate the effective-medium quantities resulting from the adjoint (or closure) problem solution and, second, to a validation through comparison between DNS of the pore-scale flow problem and predictions from the upscaled equation in a model two-dimensional configuration. Conclusions are drawn in § 5.
2. Pore-scale flow model
The development starts by considering the incompressible, Newtonian, isothermal and creeping flow of two immiscible fluid phases 
$\beta$ and 
$\gamma$ (indifferently denoted the 
$\alpha$-phase in the following) through a rigid and homogeneous porous medium, the solid skeleton of which is the 
$\sigma$-phase. No mass transport is assumed to take place at the fluid–fluid, 
$\mathscr {A}_{\beta \gamma }$, and solid–fluid, 
$\mathscr {A}_{\alpha \sigma }$, interfaces. Both fluid phases are assumed to occupy a connected region of the pore space, and the analysis is focused on a periodic unit cell, which is assumed to be representative of the process and implies that both 
$\mathscr {A}_{\beta \gamma }$ and 
$\mathscr {A}_{\alpha \sigma }$ are periodic. The flow is driven by volume forces per unit mass, 
$\boldsymbol {b}_\alpha$, in each phase and by macroscopic pressure gradients applied at the outer boundaries of the unit cell, 
$\boldsymbol {\nabla }\langle p_\alpha \rangle ^\alpha$, all being considered constant in the remainder of this work. The pore-scale velocities can therefore be considered as periodic fields. The flow problem, in dimension 
$N$ of space, can be written in a unit cell as (
$\alpha =\beta,\gamma$)
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_\alpha =0, \quad \text{in } \mathscr{V}_\alpha, \end{gather}
\begin{gather}\textbf{0}= \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{p_\alpha}+\rho_\alpha \boldsymbol{b}_\alpha , \quad \text{in } \mathscr{V}_\alpha, \end{gather}
\begin{gather}{\boldsymbol{v}}_\beta = {\boldsymbol{v}}_\gamma, \quad \text{at } \mathscr{A}_{\beta\gamma}, \end{gather}
\begin{gather}{\boldsymbol{n}}_{\beta \gamma} \boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{p_\beta}-{\boldsymbol{\mathsf{T}}}_{p_\gamma}) = 2H \gamma {\boldsymbol{n}}_{\beta \gamma}+\boldsymbol{\nabla}_s \gamma, \quad \text{at } \mathscr{A}_{\beta\gamma}, \end{gather}
\begin{gather}{\boldsymbol{v}}_\alpha =\textbf{0}, \quad \text{at } \mathscr{A}_{\alpha \sigma}, \end{gather}
\begin{gather}\boldsymbol{v}_\alpha \big|_{\mathscr{S}_{\alpha i}^-} =\boldsymbol{v}_\alpha \big|_{\mathscr{S}_{\alpha i}^+}, \quad i=1,\dots, N, \end{gather}
\begin{gather}(\boldsymbol{n}_\alpha\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{p_\alpha}) _{\mathscr{S}_{\alpha i}^-}={-}(\boldsymbol{n}_\alpha\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{p_\alpha}) _{\mathscr{S}_{\alpha i}^+} -\boldsymbol{l}_i\boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha \boldsymbol{n}_\alpha\big|_{\mathscr{S}_{\alpha i}^+}, \quad i=1,\dots, N, \end{gather}
\begin{gather}p_\alpha=0, \quad \text{at } \boldsymbol{r}_\alpha^0. \end{gather} In the above equations, 
$\mathscr {V}_\alpha$ is the space (of measure 
$V_\alpha$) occupied by the 
$\alpha$-phase within the unit cell of volume 
$V$. In addition, 
$\rho _\alpha$, 
${\boldsymbol {v}}_\alpha$ and 
${\boldsymbol{\mathsf{T}}}_{p_\alpha } = -p_\alpha {\boldsymbol{\mathsf{I}}}+ \mu _\alpha ( \boldsymbol {\nabla } {\boldsymbol {v}}_\alpha +\boldsymbol {\nabla }{\boldsymbol {v}}_\alpha ^{\rm T})$, respectively, represent the fluid density, velocity vector and total stress tensor in each phase, with 
$p_\alpha$ and 
$\mu _\alpha$, respectively, denoting the pore-scale 
$\alpha$-phase pressure and dynamic (constant) viscosity, whereas 
${\boldsymbol{\mathsf{I}}}$ is the identity tensor. In (2.1d), the stress jump is due to capillary effects, including a possible surface gradient, 
$\boldsymbol {\nabla }_s\gamma$, of the interfacial tension, 
$\gamma$, along 
$\mathscr {A}_{\beta \gamma }$ at which the mean curvature and unit normal vector, directed from the 
$\beta$-phase towards the 
$\gamma$-phase, are respectively denoted by 
$H$ and 
${\boldsymbol {n}}_{\beta \gamma }$. Note that (2.1h) is required for the problem to be well posed. Moreover, the flow is supposed to be free of acceleration, yielding a steady version of the momentum equation (2.1b). Nevertheless, flow unsteadiness can arise due to the displacement of the fluid–fluid interfaces, with or without any change of the fluids’ volume fraction, which, in the former case, justifies the unsteady character of the macroscopic mass balance equation (1.1b) and, in the latter, yields a flow that is macroscopically steady.
In (2.1f) and (2.1g), 
$\mathscr {S}_{\alpha i}^-$ and 
$\mathscr {S}_{\alpha i}^+$ represent the unit cell outer boundaries, so that 
$\boldsymbol {r}_\alpha \big |_{\mathscr {S}_{\alpha i}^+}=\boldsymbol {r}_\alpha \big |_{\mathscr {S}_{\alpha i}^-}+\boldsymbol {l}_i$, where 
$\boldsymbol {l}_i$ is the unit cell lattice vector in the 
$i{{\rm th}}$ direction (
$\alpha =\beta,\gamma$). In the rest of this work, 
$\boldsymbol {r}_{\alpha }$ is used to locate a point within the 
$\alpha$-phase with respect to a fixed coordinate system. This vector can be decomposed as 
${\boldsymbol {r}}_\alpha = {\boldsymbol {x}}_\alpha +{\boldsymbol {z}}_\alpha$, with 
${\boldsymbol {x}}_\alpha$ locating the 
$\alpha$-phase barycentre and 
${\boldsymbol {z}}_\alpha$ locating the same point as 
${\boldsymbol {r}}_\alpha$, albeit with respect to 
${\boldsymbol {x}}_\alpha$. Furthermore, the intrinsic average, 
$\langle \psi _\alpha \rangle ^\alpha$, of a pore-scale quantity, 
$\psi _\alpha$, locates the resulting average at the phase barycentre, i.e.
\begin{equation} \frac{1}{V_\alpha}\int_{\mathscr{V}_\alpha} \psi_\alpha\,{\rm d}V= \langle \psi_\alpha\rangle^\alpha\big|_{\boldsymbol{x}_\alpha}, \quad \alpha=\beta,\gamma. \end{equation} Note that 
${\boldsymbol {x}}_\alpha = \langle {\boldsymbol {r}}_\alpha \rangle ^\alpha$ and, consequently, 
$\langle \boldsymbol {z}_\alpha \rangle ^\alpha =\boldsymbol {0}$. In addition, the superficial average, 
$\langle \psi _\alpha \rangle _\alpha$, of 
$\psi _\alpha$ is defined as 
$\langle \psi _\alpha \rangle _\alpha = \varepsilon _\alpha \langle \psi _\alpha \rangle ^\alpha$, where 
$\varepsilon _\alpha =V_\alpha /V$ is the 
$\alpha$-phase volume fraction. Following Gray (Reference Gray1975), 
$\psi _\alpha$ can be decomposed in terms of its intrinsic average and spatial deviations, 
$\tilde {\psi }_\alpha$, according to
\begin{equation} \psi_\alpha\big|_{\boldsymbol{r}_\alpha} =\tilde{\psi}_\alpha\big|_{\boldsymbol{r}_\alpha} +\langle\psi_\alpha\rangle^\alpha\big|_{\boldsymbol{r}_\alpha}, \quad \alpha=\beta,\gamma. \end{equation} Note that 
$\tilde {\psi }_\alpha$ is 
$\boldsymbol {l}_i$-periodic. Performing a Taylor-series expansion about 
${\boldsymbol {x}}_\alpha$ of the above decomposition and averaging the result yields (
$\alpha =\beta,\gamma$)
\begin{align} \langle \psi_\alpha\rangle^\alpha \big|_{{\boldsymbol{x}}_\alpha} &=\big\langle \langle \psi_\alpha\rangle^\alpha\big|_{{\boldsymbol{x}}_\alpha} +{\boldsymbol{z}}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha \big|_{{\boldsymbol{x}}_\alpha} +\tfrac{1}{2} {\boldsymbol{z}}_\alpha{\boldsymbol{z}}_\alpha \boldsymbol{:} \boldsymbol{\nabla}\boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha \big|_{{\boldsymbol{x}}_\alpha} +\cdots \big\rangle^\alpha \big|_{\boldsymbol{x}_\alpha}\nonumber\\ &\quad +\langle\tilde{\psi}_\alpha\rangle^\alpha\big|_{\boldsymbol{x}_\alpha}. \end{align}
Since 
$\langle \psi _\alpha \rangle ^\alpha \big |_{{\boldsymbol {x}}_\alpha }$ and its successive gradients are constant within the unit cell, they can be taken out of the first average on the right-hand side of this expression, and making use of the order-of-magnitude estimate 
$\frac {1}{2}{\boldsymbol {z}}_\alpha {\boldsymbol {z}}_\alpha \boldsymbol {:} \boldsymbol {\nabla \nabla } \langle \psi _\alpha \rangle ^\alpha \big |_{{\boldsymbol {x}}_\alpha }= \boldsymbol {O}(r_0^2/L^2\langle \psi _\alpha \rangle ^\alpha )$, where 
$r_0$ and 
$L$ are the characteristic sizes of the averaging and macroscopic domains, respectively, it can be concluded from the above expression that 
$\langle \tilde {\psi }_\alpha \rangle ^\alpha \simeq 0$, as long as 
$r_0^2\ll L^2$. This average constraint serves the purpose of bounding the fields of 
$\tilde {\psi }$ and is necessary to make the problem well-posed, as shown below.
The decomposition defined in (2.3) can be applied to the pore-scale pressure in both phases, so that (2.1b), (2.1d), (2.1g) and (2.1h), respectively, can be replaced by (
$\alpha =\beta,\gamma$)
\begin{gather} \textbf{0}= \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha}- \boldsymbol{\nabla} \mathscr{P}_\alpha , \quad \text{in } \mathscr{V}_\alpha, \end{gather}
\begin{gather}{\boldsymbol{n}}_{\beta \gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\beta} ={\boldsymbol{n}}_{\beta \gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} +{\boldsymbol{n}}_{\beta \gamma} ( \langle p_\beta \rangle^\beta - \langle p_\gamma \rangle^\gamma )_{\boldsymbol{r}_{\beta \gamma}} + 2H \gamma {\boldsymbol{n}}_{\beta \gamma}+\boldsymbol{\nabla}_s \gamma, \quad \text{at } \mathscr{A}_{\beta\gamma}, \end{gather}
\begin{gather}(\boldsymbol{n}_\alpha\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha})_{\mathscr{S}_{\alpha i}^-}={-}(\boldsymbol{n}_\alpha\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha})_{\mathscr{S}_{\alpha i}^+}, \quad i=1,\dots, N, \end{gather}
\begin{gather}\langle \tilde{p}_\alpha\rangle^\alpha=0. \end{gather}
Here, 
${\boldsymbol{\mathsf{T}}}_{\tilde {p}_\alpha } = -\tilde {p}_\alpha {\boldsymbol{\mathsf{I}}}+ \mu _\alpha ( \boldsymbol {\nabla } {\boldsymbol {v}}_\alpha +\boldsymbol {\nabla } {\boldsymbol {v}}_\alpha ^{\rm T} )$ and 
$\boldsymbol {\nabla } \mathscr {P}_\alpha = \boldsymbol {\nabla } \langle p_\alpha \rangle ^\alpha -\rho _\alpha \boldsymbol {b}_\alpha$, yielding the formulation used in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) in which 
$\boldsymbol {b}_\alpha$ was taken as the gravitational acceleration. In (2.5b), 
$\boldsymbol {r}_{\beta \gamma }\equiv \boldsymbol {r}_{\beta }=\boldsymbol {r}_{\gamma }$ is used to locate points at 
$\mathscr {A}_{\beta \gamma }$. Note that the average constraint expressed in (2.5d) replaces (2.1h). In the following, 
$\boldsymbol {\nabla } \mathscr {P}_\alpha$ is treated as a constant, and a formalism similar to that reported in the above-cited reference is employed to derive an expression for the average pressure difference between the two fluid phases.
3. Derivation of the average pressure difference
At this point, it is convenient to introduce a new pair of closure variables in each phase, 
$f_\alpha$ and 
$\boldsymbol {f}_\alpha$, which solve the following adjoint (or closure) boundary-value problem in a periodic unit cell (
$\alpha =\beta,\gamma$):
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{f}}_\beta =\frac{1}{V_\beta}, \quad \text{in } \mathscr{V}_\beta, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{f}}_\gamma ={-}\frac{1}{V_\gamma}, \quad \text{in } \mathscr{V}_\gamma, \end{gather}
\begin{gather}\textbf{0}= \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\alpha}, \quad \text{in }\mathscr{V}_\alpha, \end{gather}
\begin{gather}{\boldsymbol{f}}_\beta = {\boldsymbol{f}}_\gamma, \quad {\boldsymbol{n}}_{\beta \gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\beta} ={\boldsymbol{n}}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\gamma}, \quad \text{at } \mathscr{A}_{\beta\gamma}, \end{gather}
\begin{gather}{\boldsymbol{f}}_\alpha =\textbf{0}, \quad \text{at } \mathscr{A}_{\alpha \sigma}, \end{gather}
\begin{gather}\boldsymbol{f}_\alpha\big|_{\mathscr{S}_{\alpha i}^-}= \boldsymbol{f}_\alpha\big|_{\mathscr{S}_{\alpha i}^+}, \quad (\boldsymbol{n}_\alpha\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{f_\alpha}) _{\mathscr{S}_{\alpha i}^-}={-}(\boldsymbol{n}_\alpha\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{f_\alpha}) _{\mathscr{S}_{\alpha i}^+}, \quad i=1,\dots, N, \end{gather}
\begin{gather}f_\alpha=0, \quad \text{at } \boldsymbol{r}_\alpha^ 0. \end{gather} In these equations, 
${\boldsymbol{\mathsf{T}}}_{f_\alpha }=-f_\alpha {\boldsymbol{\mathsf{I}}} +\mu _\alpha (\boldsymbol {\nabla }\boldsymbol {f}_\alpha +\boldsymbol {\nabla } \boldsymbol {f}_\alpha ^{\rm T})$ is a stress-like second-order tensor, whereas, in (3.1i), 
$\boldsymbol {r}_\alpha ^0$ locates an arbitrary point in the 
$\alpha$-phase (
$\alpha =\beta,\gamma$). This new closure problem can be combined with the flow problem in a unit cell (namely (2.1a), (2.1c), (2.5), (2.1e) and (2.1f)) by means of a Green's formula. Taking two arbitrary scalar fields, 
$a$ and 
$q$, and two arbitrary vector fields, 
$\boldsymbol {a}$ and 
$\boldsymbol {q}$, all of them defined in the 
$\alpha$-phase and having sufficient regularity, this Green's formula can be written as (see the derivation in Appendix A of Sánchez-Vargas, Valdés-Parada & Lasseux (Reference Sánchez-Vargas, Valdés-Parada and Lasseux2022))
\begin{align} & \int_{\mathscr{V}_\alpha} [ {\boldsymbol{a}}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_q)- (\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_a)\boldsymbol{\cdot} {\boldsymbol{q}} -q\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{a}} + a \boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{q}}]\,{\rm d}V \nonumber\\ & \quad = \int_{\mathscr{A}_\alpha} [ {\boldsymbol{a}}\boldsymbol{\cdot} ({\boldsymbol{n}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_q) - {\boldsymbol{n}}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_a\boldsymbol{\cdot} {\boldsymbol{q}} ]\,{\rm d}A. \end{align} Here 
$\mathscr {A}_\alpha$ represents the enclosing surfaces of 
$\mathscr {V}_\alpha$, whereas 
${\boldsymbol{\mathsf{T}}}_a\equiv -a{\boldsymbol{\mathsf{I}}}+\mu _\alpha (\boldsymbol {\nabla }\boldsymbol {a} +\boldsymbol {\nabla }\boldsymbol {a}^{\rm T})$ and 
${\boldsymbol{\mathsf{T}}}_q\equiv -q{\boldsymbol{\mathsf{I}}} +\mu _\alpha (\boldsymbol {\nabla }\boldsymbol {q} +\boldsymbol {\nabla }\boldsymbol {q}^{\rm T})$ (
$\alpha =\beta,\gamma$). Fixing 
$a\equiv \tilde {p}_\beta$, 
$\boldsymbol {a}\equiv \boldsymbol {v}_\beta$, 
$q\equiv f_\beta$, 
$\boldsymbol {q}\equiv \boldsymbol {f}_\beta$ and 
$\mu _\alpha \equiv \mu _\beta$ in the above equation, and substituting the corresponding differential equations and boundary conditions, together with the constraint (2.5d), leads to the following expression:
\begin{align} -\int_{\mathscr{V}_\beta}\boldsymbol{f}_\beta \,{\rm d}V \boldsymbol{\cdot}\boldsymbol{\nabla} \mathscr{P}_\beta &=\int_{\mathscr{A}_{\beta\gamma}} [ {\boldsymbol{v}}_\gamma \boldsymbol{\cdot} ( {\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\gamma} ) - ({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}) \boldsymbol{\cdot}{\boldsymbol{f}}_\gamma ]\,{\rm d}A \nonumber\\ &\quad -\int_{\mathscr{A}_{\beta\gamma}} ( \langle p_\beta\rangle^\beta -\langle p_\gamma\rangle^\gamma)_{\boldsymbol{r}_{\beta\gamma}}{\boldsymbol{n}}_{\beta \gamma}\boldsymbol{\cdot}{\boldsymbol{f}}_\beta\,{\rm d}A\nonumber\\ &\quad -\int_{\mathscr{A}_{\beta\gamma}} (2 H \gamma {\boldsymbol{n}}_{\beta \gamma}+\boldsymbol{\nabla}_s \gamma ) \boldsymbol{\cdot}{\boldsymbol{f}}_\beta \,{\rm d}A. \end{align} Note that 
$\boldsymbol {\nabla }\mathscr {P}_\beta$ is taken out of the integral on the left-hand side, as it is assumed to be constant. However, this is not the case for the average pressure difference inside the second area integral on the right-hand side of this equation, as both average pressures are considered at points 
$\boldsymbol {r}_{\beta \gamma }$ at 
$\mathscr {A}_{\beta \gamma }$.
To progress towards an expression for the average pressure difference, a Taylor-series expansion of 
$\langle p_\alpha \rangle ^\alpha|_{\boldsymbol {r}_{\beta \gamma }}$ can be performed about 
${\boldsymbol {x}}_\alpha$. Maintaining the first two terms of this expansion (the remaining ones are indeed zero if 
$\boldsymbol {\nabla } \langle p_\alpha \rangle ^\alpha$ is constant), and recalling that 
$\boldsymbol {\nabla } \mathscr {P}_\beta = \boldsymbol {\nabla } \langle p_\beta \rangle ^\beta -\rho _\beta \boldsymbol {b}_\beta$, the above equation takes the following form:
\begin{align} &-\int_{\mathscr{V}_\beta}{\boldsymbol{f}}_\beta
\,{\rm d}V \boldsymbol{\cdot} (\boldsymbol{\nabla} \langle
p_\beta\rangle^\beta-\rho_\beta\boldsymbol{b}_\beta)
\nonumber\\ &\quad =
\int_{\mathscr{A}_{\beta\gamma}} [
{\boldsymbol{v}}_\gamma \boldsymbol{\cdot}
({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot}
{\boldsymbol{\mathsf{T}}}_{f_\gamma}) -
({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot}
{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma})\boldsymbol{\cdot}
{\boldsymbol{f}}_\gamma ]\,{\rm d}A \nonumber\\ &\qquad - (\langle
p_\beta\rangle^\beta\big|_{\textbf{x}_\beta} -\langle
p_\gamma\rangle^\gamma\big|_{{\boldsymbol{x}}_\gamma})
\int_{\mathscr{A}_{\beta\gamma}} {\boldsymbol{n}}_{\beta
\gamma} \boldsymbol{\cdot}{\boldsymbol{f}}_\beta\,{\rm d}A
-{\int}_{\mathscr{A}_{\beta\gamma}} {\boldsymbol{n}}_{\beta
\gamma}\boldsymbol{\cdot} {\boldsymbol{f}}_\beta
{\boldsymbol{z}}_\beta\,{\rm d}A \boldsymbol{\cdot}
\boldsymbol{\nabla} \langle p_\beta\rangle^\beta
\nonumber\\ &\qquad +
\int_{\mathscr{A}_{\beta\gamma}}{\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot}{\boldsymbol{f}}_\beta{\boldsymbol{z}}_\gamma\,{\rm
d}A \boldsymbol{\cdot}\boldsymbol{\nabla} \langle
p_\gamma\rangle^\gamma -{\int_{\mathscr{A}_{\beta\gamma}}
(}2 H \gamma {\boldsymbol{n}}_{\beta
\gamma}+\boldsymbol{\nabla}_s \gamma )
\boldsymbol{\cdot}{\boldsymbol{f}}_\beta\,{\rm d}A.
\end{align}
Note that in the area integrals on the right-hand side of this expression, 
${\boldsymbol {z}}_\alpha = {\boldsymbol {r}}_{\beta \gamma } - {\boldsymbol {x}}_\alpha$ (
$\alpha =\beta,\gamma$).
To simplify the above equation, it is pertinent to note that 
$\int _{\mathscr {A}_{\beta \gamma }}{\boldsymbol {n}}_{\beta \gamma }\boldsymbol {\cdot } {\boldsymbol {f}}_\beta \,{\rm d}A=$ 
$\int _{\mathscr {V}_\beta } \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {f}}_\beta \,{\rm d}V =1$. In addition, the following nomenclature is proposed for the effective-medium coefficients (
$\alpha =\beta,\gamma$ and 
$\alpha \ne \kappa$)
\begin{gather} \boldsymbol{\varphi}_\alpha = \int_{\mathscr{V}_\alpha} \boldsymbol{f}_\alpha\,{\rm d}V, \quad \boldsymbol{\psi}_\alpha=\boldsymbol{\varphi}_\alpha -\int_{\mathscr{A}_{\beta\gamma}} {\boldsymbol{n}}_{\alpha \kappa}\boldsymbol{\cdot} {\boldsymbol{f}}_\alpha \boldsymbol{z}_\alpha\,{\rm d}A, \end{gather}
\begin{gather} s_{\beta\gamma}= \int_{\mathscr{A}_{\beta \gamma}} (2H\gamma {\boldsymbol{n}}_{\beta \gamma}+\boldsymbol{\nabla}_s \gamma) \boldsymbol{\cdot} {\boldsymbol{f}}_{\beta}\,{\rm d}A,\end{gather}and for the average pressure difference
\begin{equation} \Delta\mathcal{P} = \langle p_\gamma\rangle^\gamma\big|_{\boldsymbol{x}_\gamma} -\langle p_\beta\rangle^\beta\big|_{\boldsymbol{x}_\beta}. \end{equation}After rearranging, (3.4) can be written in the following more compact form:
\begin{align} \Delta\mathcal{P} &={-}\int_{\mathscr{A}_{\beta\gamma}} [{\boldsymbol{v}}_\gamma \boldsymbol{\cdot} ({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\gamma}) - {\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}\boldsymbol{\cdot} {\boldsymbol{f}}_\gamma]\,{\rm d}A -\boldsymbol{\psi}_\beta \boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\beta\rangle^\beta \nonumber\\ &\quad + \rho_\beta {\boldsymbol{b}}_\beta\boldsymbol{\cdot}\boldsymbol{\varphi}_\beta -(\boldsymbol{\psi}_\gamma -\boldsymbol{\varphi}_\gamma) \boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma +s_{\beta\gamma}. \end{align}
The first area integral in (3.7) can be expressed in an alternative form by considering again Green's formula given in (3.2) with 
$a\equiv \tilde {p}_\gamma$, 
${\boldsymbol {a}}\equiv {\boldsymbol {v}}_\gamma$, 
$q\equiv f_\gamma$, 
${\boldsymbol {q}}\equiv {\boldsymbol {f}}_\gamma$ and 
$\mu _\alpha \equiv \mu _\gamma$. After substitution of the corresponding differential equations and boundary conditions, the resulting expression is
\begin{equation} \int_{\mathscr{A}_{\beta\gamma}} [\textbf{v}_\gamma\boldsymbol{\cdot} ({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{f_\gamma}) -({\boldsymbol{n}}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma})\boldsymbol{\cdot} {\boldsymbol{f}}_\gamma ]\,{\rm d}A = \boldsymbol{\varphi}_\gamma \boldsymbol{\cdot} \boldsymbol{\nabla} \mathscr{P}_\gamma. \end{equation}Substituting this result into (3.7) yields the following expression for the average pressure difference:
\begin{equation} \Delta\mathcal{P}={-}\boldsymbol{\psi}_\beta \boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\beta\rangle^\beta +\rho_\beta \boldsymbol{b}_\beta\boldsymbol{\cdot} \boldsymbol{\varphi}_\beta -\boldsymbol{\psi}_\gamma\boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma +\rho_\gamma \boldsymbol{b}_\gamma\boldsymbol{\cdot}\boldsymbol{\varphi}_\gamma +s_{\beta \gamma}. \end{equation}An alternative equivalent form can be derived by combining the adjoint problem defined in (3.1) with the flow problem given in (2.1) by means of Green's formula (3.2). Employing the same procedure as above results in the following expression:
\begin{align} \Delta\mathcal{P} &= \sum_{i=1}^N\left[\int_{\mathscr{S}_{\beta i}^+} \boldsymbol{n}_\beta\boldsymbol{\cdot}\boldsymbol{f}_\beta \,{\rm d}A\, \boldsymbol{l}_i\boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\beta \rangle^\beta + \int_{\mathscr{S}_{\gamma i}^+} \boldsymbol{n}_\gamma\boldsymbol{\cdot}\boldsymbol{f}_\gamma \,{\rm d}A\, \boldsymbol{l}_i\boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\gamma \rangle^\gamma\right] \nonumber\\ &\quad + \rho_\beta \boldsymbol{b}_\beta\boldsymbol{\cdot} \boldsymbol{\varphi}_\beta +\rho_\gamma\boldsymbol{b}_\gamma \boldsymbol{\cdot} \boldsymbol{\varphi}_\gamma +s_{\beta\gamma}. \end{align}
Using the identity 
$\int _{\mathscr {V}_\alpha }\boldsymbol {\nabla } \boldsymbol {\cdot }(\,\boldsymbol {f}_\alpha \boldsymbol {z}_\alpha )\,{\rm d}V = \int _{\mathscr {A}_\alpha }\boldsymbol {n}_\alpha \boldsymbol {\cdot } \boldsymbol {f}_\alpha \boldsymbol {z}_\alpha \,{\rm d}A = \int _{\mathscr {V}_\alpha }\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {f}_\alpha \boldsymbol {z}_\alpha \,{\rm d}V + \int _{\mathscr {V}_\alpha }\boldsymbol {f}_\alpha \,{\rm d}V$ (note that 
$\boldsymbol {\nabla }\boldsymbol {z}_\alpha ={\boldsymbol{\mathsf{I}}}$ is employed here), taking into account (3.1a,b) and boundary conditions for 
${\boldsymbol {f}}_\alpha$, it can be readily shown that 
$\sum _{i=1}^N\int _{\mathscr {S}_{\alpha i}^+} \boldsymbol {n}_\alpha \boldsymbol {\cdot }\boldsymbol {f}_\alpha \,{\rm d}A\, \boldsymbol {l}_i=-\boldsymbol {\psi }_\alpha$ (
$\alpha =\beta,\gamma$), proving that (3.9b) coincides with (3.9a). Note that the average pressure difference can also be expressed at the centroid, 
$\boldsymbol {x}$, of the unit cell by using the relationship 
$\langle p_\alpha \rangle ^\alpha \big |_{\boldsymbol {x}}= \langle p_\alpha \rangle ^\alpha \big |_{\boldsymbol {x}_\alpha } -\langle \boldsymbol {y}_\alpha \rangle ^\alpha \boldsymbol {\cdot } \boldsymbol {\nabla }\langle p_\alpha \rangle ^\alpha$, with the definition 
$\boldsymbol {r}_\alpha =\boldsymbol {x}+\boldsymbol {y}_\alpha$. Indeed, equality of the two pressure gradients in both phases is a compatibility requirement with the assumption of a periodic unit cell representative of the process (see details in Appendix B in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022)). The effective-medium coefficient associated with the pressure gradient in that case can be shown to be 
$-\boldsymbol {\varphi }_\beta -\boldsymbol {\varphi }_\gamma +\boldsymbol {x}_\gamma -\boldsymbol {x}_\beta$.
The expressions for 
$\Delta \mathcal {P}$ given in (3.9), which are the important result of this work, only rely on the assumptions made for the pore-scale flow problem together with scale separation and the possibility of a local flow description in a representative periodic unit cell. These assumptions may appear to be overly severe and were adopted with the intention of deriving a closure problem for the computation of the effective-medium quantities involved in (3.9). Nevertheless, the upscaled expressions may be used even in situations when these assumptions are not completely met. Note that the derivation does not require 
$\mathscr {A}_{\beta \gamma }$ to be stationary (i.e. motionless fluid–fluid interfaces), or a time-independent saturation, as explained in § 2. The first four terms on the right-hand side of (3.9a) represent the influence of the hydrodynamic forcing in each phase, whereas the last term represents the influence of interfacial tension. The latter results from a complex interfacial average of the Laplace effects (i.e. 
$2H\gamma$) and interfacial tension gradient effects at the pore scale, and involves 
$\boldsymbol {f}_\beta \ (=\boldsymbol {f}_\gamma )$ at 
$\mathscr {A}_{\beta \gamma }$ as a weighting function.
Obviously, 
$\Delta \mathcal {P}$ depends not only on the wetting-phase saturation, but also on the configuration of 
$\mathscr {A}_{\beta \gamma }$ as well as on the dynamic effects, as anticipated, for instance, by Hassanizadeh & Gray (Reference Hassanizadeh and Gray1993) and later by Gray & Miller (Reference Gray and Miller2011), but yet in an unclosed form. This is clearly different from the Laplace equation typically found in the literature (Whitaker Reference Whitaker1994). This result is only recovered in the absence of any macroscopic forcing. Indeed, in that case, 
$2H$ must be constant, and, assuming that 
$\gamma$ is also constant, (3.9a) reduces to the classical macroscopic Laplace relation, i.e. 
$\Delta \mathcal {P}\equiv \mathcal {P}_{cL}=2H\gamma$. The vectors 
$\boldsymbol {\psi }_\alpha$ and 
$\boldsymbol {\varphi }_\alpha$ can be conceived as effective characteristic lengths for, respectively, 
$\boldsymbol {\nabla }\langle p_\alpha \rangle ^\alpha$ and 
$\rho _\alpha {\boldsymbol {b}}_\alpha$.
The effective-medium coefficients, 
$\boldsymbol {\psi }_\alpha$, 
$\boldsymbol {\varphi }_\alpha$ and 
$s_{\beta \gamma }$, can be predicted from the solution of the associated closure problem defined in (3.1). The equation for 
$\Delta \mathcal {P}$ is compatible with an expression proposed in the literature (see (43) in Hassanizadeh & Gray (Reference Hassanizadeh and Gray1993)). It has the advantage of being closed and it allows identifying and evaluating the role of each source. An analysis of the prediction from this upscaled equation is proposed in the next section.
4. Results
In this section, numerical simulation results are reported for a two-dimensional (
$N=2$) two-phase flow configuration schematically represented in figure 1 and already employed by Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) to investigate the macroscopic momentum equations. Results are presented for 
$\boldsymbol {b}_\alpha =\boldsymbol {0}$, 
$\boldsymbol {\nabla }_s\gamma =\boldsymbol {0}$ and a prescribed pressure gradient in each phase given by 
$\partial \langle p_\alpha \rangle ^\alpha /\partial x=-h$ and 
$\partial \langle p_\alpha \rangle ^\alpha /\partial y=0$ (
$\alpha =\beta,\gamma$). Moreover, the flow is considered to be fully steady, i.e. with 
$\mathscr {A}_{\beta \gamma }$ at its stationary position for the sake of brevity in presentation.
Two-dimensional two-phase flow configuration considered for the numerical simulations, made of a square pattern of parallel cylinders of circular cross-section (a) and corresponding periodic unit cell (b). The flow results from a constant horizontal macroscopic pressure gradient in each phase. The 
$\beta$-phase flows as a wetting layer attached to each cylinder row, whereas the non-wetting 
$\gamma$-phase flows in the intermediate layers.
4.1. Flow and closure problem solutions
The first step in the analysis is the solution of the pore-scale flow problem, i.e. DNS and the computation of the average pressure difference. The boundary integral element method (BIEM) employed for DNS and results on the configuration under concern are detailed in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022) (see Appendix C for the BIEM and § 6 for the results). Accordingly, the same dimensionless quantities as in this reference are employed, i.e. 
$\boldsymbol {r}_\alpha ^*=\boldsymbol {r}_\alpha /\ell$, 
$\boldsymbol {\nabla }^*=\ell \boldsymbol {\nabla }$, 
$\boldsymbol {v}_\alpha ^*=\mu _\gamma \boldsymbol {v}_\alpha /(h\ell ^2)$ and 
$p_\alpha ^*=\mu _\gamma p_\alpha /(\mu _\alpha h\ell )$ (
$\alpha =\beta,\gamma$), where 
$\ell$ is the 
$x$ and 
$y$ unit cell size (see figure 1). The solution depends on the porosity, 
$\varepsilon =(V_\beta +V_\gamma )/V$, wetting-phase saturation, 
$S_\beta =\varepsilon _\beta /\varepsilon = V_\beta /(V_\beta +V_\gamma )$, capillary number, 
$Ca=h\ell ^2/\gamma$, and viscosity ratio, 
$\mu ^*=\mu _\beta /\mu _\gamma$. Once the steady shape of 
$\mathscr {A}_{\beta \gamma }$ is determined, the pressure fields are computed (see (C10) in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022)), together with their intrinsic average, 
$\langle p_\alpha \rangle ^\alpha \big |_{\boldsymbol {x}_\alpha }$, and their difference, yielding the DNS result for this quantity.
In the second step, the closure problem is solved taking the shape of 
$\mathscr {A}_{\beta \gamma }$ determined from DNS in the unit cell. Using the following dimensionless changes of variables, 
$d_\alpha ^*=({\ell ^3}/{\mu _\alpha })f_\alpha$, 
$\boldsymbol {d}_\beta ^*=\ell ^2(\,\boldsymbol {f}_\beta - {\boldsymbol {r}_\beta }/{N V_\beta })=\boldsymbol {f}_\beta ^*- {\boldsymbol {r}_\beta ^*}/{N\varepsilon _\beta }$ and 
$\boldsymbol {d}_\gamma ^*=\ell ^2(\,\boldsymbol {f}_\gamma + {\boldsymbol {r}_\gamma }/{NV_\gamma })= \boldsymbol {f}_\gamma ^*+{\boldsymbol {r}_\gamma ^*}/{N\varepsilon _\gamma }$, and noticing that 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {r}_\alpha =N$ and 
$\boldsymbol {\nabla }\boldsymbol {r}_\alpha ={\boldsymbol{\mathsf{I}}}$, 
$d_\alpha ^*$ and 
$\boldsymbol {d}_\alpha ^*$ satisfy a homogeneous incompressible Stokes-like problem along with the corresponding boundary conditions. Its solution is obtained with the same BIEM as for DNS. The solution on 
$\boldsymbol {d}_\alpha ^*$ (and hence on 
$\boldsymbol {f}_\alpha ^*$) at the boundaries is then used to predict the dimensionless average pressure difference, 
$\Delta \mathcal {P}^*= \langle p_\gamma ^*\rangle ^\gamma \big |_{\boldsymbol {x}^*_\gamma }-\mu ^* \langle p_\beta ^*\rangle ^\beta \big |_{\boldsymbol {x}^*_\beta }$, from the upscaled equation given by (3.9b). Under the conditions employed here, it reduces to
\begin{equation} \Delta\mathcal{P}^* ={-}\int_{\mathscr{S}_{\gamma 1}^+}f_{\gamma x}^*\,{\rm d}A^* -\int_{\mathscr{S}_{\beta 1}^+}f_{\beta x}^*\,{\rm d}A^* +s_{\beta\gamma}^*,\end{equation}
where 
$s_{\beta \gamma }^*=({1}/{Ca})\int _{\mathscr {A}_{\beta \gamma }}2H^* {\boldsymbol {n}}_{\beta \gamma }\boldsymbol {\cdot }\boldsymbol {f}_\beta ^*\,{\rm d}A^*$, 
$f_{\alpha x}^*$ is the 
$x$ component of 
$\boldsymbol {f}_\alpha ^*$ (
$\alpha =\beta,\gamma$) and 
$H^*=\ell H$.
Numerical simulations are performed for a range of 
$S_\beta$ with 
$\varepsilon =0.8$, 
$Ca=0.1$ and 
$10$, and 
$\mu ^*=0.1$ and 
$10$. The numerical results obtained with the BIEM were verified to match those obtained using the finite-element method with Comsol Multiphysics 6.0.
4.2. Validation with numerical simulations
Before focusing on the predictions of 
$\Delta \mathcal {P}^*$, it is pertinent to examine the dimensionless pressure and closure variable pore-scale fields, the boundary values of the latter providing the effective-medium quantities involved in (4.1). An illustration of 
$p_\alpha ^*$ as well as of the 
$x$ and 
$y$ components of 
$\boldsymbol {f}_\beta ^*$ and 
$\boldsymbol {f}_\gamma ^*$ are, respectively, represented in figures 2(a), 2(b) and 2(c) for 
$\varepsilon =0.8$, 
$Ca=10$, 
$\mu ^*=10$ and 
$S_\beta =0.5$. The expected symmetry of 
$p_\alpha ^*$ and 
$f_{\alpha x}^*$ and antisymmetry of 
$f_{\alpha y}^*$ (
$\alpha =\beta,\gamma$) with respect to the 
$x$-axis is clearly noticeable. The pressure jump due to surface tension and normal viscous stress contrast, and continuity of 
$\boldsymbol {f}_\alpha ^*$, are also evident. Note that, due to the non-symmetry of 
$\mathscr {A}_{\beta \gamma }$ with respect to the 
$y$-axis, no 
$y$ symmetry is observed on these fields.
(a) Dimensionless pressure field resulting from DNS at steady state; (b) 
$x$ component and (c) 
$y$ component of the dimensionless closure variable 
$\boldsymbol {f}_\alpha ^*$ (
$\alpha =\beta,\gamma$). Corresponding isolines are represented in white. Black lines materialize 
$\mathscr {A}_{\beta \sigma }$ and 
$\mathscr {A}_{\beta \gamma }$. Here 
$\varepsilon =0.8$, 
$Ca=10$, 
$\mu ^*=10$ and 
$S_\beta =0.5$.
A comparison of 
$\Delta \mathcal {P}^*$ obtained from DNS and predicted from the upscaled equation is represented versus 
$S_\beta$ in figure 3. In all cases, an excellent agreement between the two approaches is obtained over the whole range of 
$S_\beta$. Indeed, the relative error, taking the DNS result as the reference, is less than 0.8 % among all cases. This confirms the pertinence of the expression for 
$\Delta \mathcal {P}$ derived here. In figure 3, the classical macroscopic Laplace capillary pressure is also reported. The dimensionless form of this quantity is given by 
$\mathcal {P}_{cL}^*=2\langle H^*\rangle _{\beta \gamma }/Ca$, where 
$\langle H^*\rangle _{\beta \gamma }=({1}/{A_{\beta \gamma }^*}) \int _{\mathscr {A}_{\beta \gamma }}H^*\,{\rm d}A^*$, and 
$A_{\beta \gamma }^*$ is the dimensionless measure of 
$\mathscr {A}_{\beta \gamma }$. From these results, it must be emphasized that 
$\mathcal {P}_{cL}^*$ cannot represent a reasonable approximation of 
$\Delta \mathcal {P}^*$, as it is always much smaller, by at least five orders of magnitude, than 
$\Delta \mathcal {P}^*$. The dominant contribution to 
$\Delta \mathcal {P}^*$ is from 
$s_{\beta \gamma }^*$, which, however, overpredicts it. Indeed, as can be observed from figure 3, the two hydrodynamic terms remain significant, except in some specific circumstances, i.e. for 
$S_\beta \gtrsim 0.6$ for 
$\mu ^*=0.1$ and 
$Ca=10$ (figure 3a) or 
$S_\beta \lesssim 0.4$ for 
$\mu ^*=0.1$ and 
$Ca=0.1$ (figure 3b).
Dimensionless average pressure difference, 
$\Delta \mathcal {P}^*$, obtained from DNS and predicted from the upscaled equation (UE), along with the interfacial term contribution, 
$s_{\beta \gamma }^*$ (see (4.1)), on the left-hand vertical axes, and classical Laplace capillary pressure, 
$\mathcal {P}_{cL}^*$ (see text), on the right-hand vertical axes, versus the wetting-phase saturation, 
$S_\beta$, for: 
$\varepsilon =0.8$; and (a) 
$\mu ^*=0.1$, 
$Ca=10$; (b) 
$\mu ^*=0.1$, 
$Ca=0.1$; (c) 
$\mu ^*=10$, 
$Ca=10$; (d) 
$\mu ^*=10$, 
$Ca=0.1$.
It was also numerically verified that the upscaled equation prediction of 
$\Delta \mathcal {P}^*$ perfectly matches the DNS results for other shapes of 
$\mathscr {A}_{\beta \gamma }$, different from the stationary one that is characterized by a zero normal interfacial velocity. Moreover, although the flow under consideration here is macroscopically steady (
$\partial S_\beta /\partial t=0$), it clearly appears that 
$\Delta \mathcal {P}$ does not reduce to the equilibrium capillary pressure remaining as the unique term if the relationship proposed by Hassanizadeh & Gray (Reference Hassanizadeh and Gray1993) is considered. This was suggested by Gray & Miller (Reference Gray and Miller2011), who stated that this relationship ‘does not account properly for disequilibrium due to the dynamics of interfacial rearrangement at constant saturation’; the upscaled equation derived here accounts for this mechanism.
5. Conclusions
In this work, an upscaled equation was derived to predict the average pressure difference, 
$\Delta \mathcal {P}$, between two immiscible fluid phases in homogeneous porous media using a simplified version of the volume-averaging method, combined with elements of the adjoint method and Green's formula. The resulting expression shows that 
$\Delta \mathcal {P}$ includes contributions from surface tension effects, and also from macroscopic pressure gradients and volume forces acting in each phase. The derivation of the expression for 
$\Delta \mathcal {P}$ is not restricted to situations in which the fluid–fluid interface is at its stationary position. The effective-medium quantities involved in this expression are predicted from the solution of a single associated closure (adjoint) problem in a representative periodic unit cell. Predictions of 
$\Delta \mathcal {P}$ from its expression derived here were validated through comparisons with DNS, showing excellent agreement. It is shown that, for the particular porous structure and phases distribution considered here, the term related to capillary effects is dominant in the prediction of the average pressure difference (especially as the wetting-phase saturation is closer to unity), albeit the influence of the macroscopic pressure gradients cannot be disregarded.
Imaging techniques may be envisaged to obtain fluid-phase distributions in a representative domain as an input for the closure problem solution and to further investigate real configurations. Alternatively, an exhaustive numerical analysis of porous skeleton topologies and fluid-phase distributions, which may be seen as representative of a real specific configuration of interest, could be performed to propose correlations that could be implemented in standard two-phase flow simulators of common use. This would avoid the dependence on the solution of closure problems and allow one to express the upscaled equation for 
$\Delta \mathcal {P}$ in terms of macroscopic parameters only.
The average pressure difference derived in this work completes the macroscopic mass and momentum equations previously reported by Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022). Finally, the analysis provided here opens the way to study more complex flow conditions that will be addressed in future works.
Funding
Part of this work was performed while both authors were visiting professors at the University of Genoa (UniGe) in collaboration with Professor A. Bottaro. Financial support from UniGe is sincerely appreciated. Data of the numerical results presented in this work are available from the authors upon request.
Declaration of interests
The authors report no conflict of interest.
