2.1 Governing equation

The problem under consideration has already been formulated by Wilmott et al. (Reference Wilmott, Breward and Chapman2018) so only brief summary details of the main governing equation are given here. The notation used here is slightly different from the study of Wilmott et al. (Reference Wilmott, Breward and Chapman2018). The reason is that the study in question considered that a given pressure difference was applied and determined the speed of advance of the droplet as part of the solution of the problem. A capillary number was therefore defined in terms of the known applied pressure rather than the a priori unknown droplet speed. Here however we adopt a simpler approach in which the droplet speed is treated as given. The capillary number then becomes $\unicode[STIX]{x1D707}V/\unicode[STIX]{x1D6FE}$ where $\unicode[STIX]{x1D707}$ is aqueous phase viscosity, $V$ is droplet speed and $\unicode[STIX]{x1D6FE}$ is the interfacial tension between the oil and aqueous phase. This corresponds to the definition used in the original formulation of Bretherton (Reference Bretherton1961), so we denote it $Ca_{B}$ (so as not to confuse it with the alternative definition of capillary number employed by Wilmott et al. (Reference Wilmott, Breward and Chapman2018)).

Solving the problem for advance of an oil droplet requires solving for the shape of the aqueous film that separates the droplet from the capillary wall: see figure 1. In the interests of simplicity we assume a two-dimensional case (rather than an axisymmetric system) with the half-thickness of the capillary being denoted $R$ . A static drop in the capillary would then have a curvature radius exactly equal to $R$ on its front and back ends. Here however the droplet is moving, so the curvature radius can be perturbed slightly. The droplet shape can in principle be affected by the ratio between the droplet viscosity and aqueous phase viscosity. However analysis by Park & Homsy (Reference Park and Homsy1984) has suggested that, provided the viscosity ratio is no larger than order $Ca_{B}^{-1/3}$ (with $Ca_{B}\ll 1$ in problems of interest), then the viscosity ratio has no bearing on the shape.

As alluded to earlier, uniformly curved capillary static regions at the ends of the drop are joined up with an aqueous thin film across a transition region. Determining the shape of the transition region then becomes part of the required solution. In what follows, unless explicitly specified otherwise, we will employ the term ‘thin film’ to refer to the thinnest part excluding the transition region, and the term ‘film’ more generically including possibly the transition region in addition to the thinnest part. It turns out moreover that the shape of the transition region near the rear of the drop (Burgess & Foster Reference Burgess and Foster1990; Wong, Radke & Morris Reference Wong, Radke and Morris1995a ,Reference Wong, Radke and Morris b ; Giavedoni & Saita Reference Giavedoni and Saita1999; Wilmott et al. Reference Wilmott, Breward and Chapman2018) is a little more complex than the shape near the front. In the interests of simplicity therefore we focus here solely on the region at the front of the drop.

We switch to a frame of reference moving with the drop. In this frame of reference, the shape of the drop (i.e. the thickness of the transition region versus distance along it) should have attained a steady state. Hence an ordinary differential equation (rather than a partial differential equation) can be solved. In what follows, equations are cast in dimensionless form. Both thickness and distance are scaled by $R$ , pressures are scaled by $\unicode[STIX]{x1D6FE}/R$ and velocities are scaled by $\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D707}$ . The governing equation for dimensionless aqueous film thickness $h$ versus dimensionless distance $Z$ along the film (which is derived in appendix A following Wilmott et al. (Reference Wilmott, Breward and Chapman2018)) is then

(2.1) $$\begin{eqnarray}h_{ZZZ}=3Ca_{B}\frac{(h-h^{\ast })}{h^{3}}-\frac{1}{\bar{\unicode[STIX]{x1D6E4}}}\left(\frac{\unicode[STIX]{x1D70E}_{c0}^{2}+2s^{\ast }\unicode[STIX]{x1D70E}_{c0}\unicode[STIX]{x1D70E}_{o0}\cosh (h/\unicode[STIX]{x1D716}_{2})+s^{\ast \,2}\unicode[STIX]{x1D70E}_{o0}^{2}}{\sinh ^{2}(h/\unicode[STIX]{x1D716}_{2})}\right)_{Z},\end{eqnarray}$$

where $Ca_{B}$ is the capillary number, $h^{\ast }$ is the uniform film thickness outside the transition region (which is a priori unknown), $\bar{\unicode[STIX]{x1D6E4}}$ is a parameter (precise mathematical definition to be given shortly) measuring the relative importance of capillary and electrostatic terms (specifically referred to electrostatics on the capillary wall) and $\unicode[STIX]{x1D716}_{2}$ (again precise mathematical definition to be given shortly) is the ratio between a Debye screening length (Wright Reference Wright2007) and channel half-thickness $R$ .

The bracketed term on the right-hand side of (2.1) represents an electro-osmotic tension (Yang, Fu & Hwang Reference Yang, Fu and Hwang2001; Waghmare & Mitra Reference Waghmare and Mitra2008). The parameter $\unicode[STIX]{x1D70E}_{c0}$ is the net charge per adsorption site on the capillary wall, $\unicode[STIX]{x1D70E}_{o0}$ is the charge per adsorption site on the oil and $s^{\ast }$ is the ratio between the density of adsorption sites on the oil and on the wall. Note that both $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ can vary in principle between $-1$ (corresponding to all of the negatively charged adsorption sites remaining unoccupied) and $+1$ (corresponding to all adsorption sites being occupied by divalent positive ions). Moreover $h/\unicode[STIX]{x1D716}_{2}$ is the ratio between film thickness and Debye screening length. In general it is necessary to consider the full $h/\unicode[STIX]{x1D716}_{2}$ dependence of the electro-osmotic tension in (2.1), which is relatively complicated. There are however some simpler limiting cases considered in appendix B. In particular, if $h/\unicode[STIX]{x1D716}_{2}\gg 1$ electrostatic effects are screened by electrical double layer effects, and so are weak regardless of the value of $\bar{\unicode[STIX]{x1D6E4}}$ . The electro-osmotic tension then reduces to an exponential decay, a case considered by Teletzke et al. (Reference Teletzke, Davis and Scriven1988).

2.2 Rescaling the governing equation

In order to solve the governing equation (2.1), a number of manipulations are required. These are standard in the Bretherton case (ignoring electrostatic effects), but in what follows we review how they are modified in the presence of electrostatic terms. It is known that (Bretherton Reference Bretherton1961) the length and thickness of the transition region scale respectively as $Ca_{B}^{1/3}$ and $Ca_{B}^{2/3}$ , scalings that apply at least nominally in the Bretherton case for which electrostatic effects are neglected. It is then convenient to define new variables $\unicode[STIX]{x1D6F9}=Z/(3Ca_{B})^{1/3}$ and $J=h/(3Ca_{B})^{2/3}$ . Equation (2.1) then becomes

(2.2) $$\begin{eqnarray}J_{\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6F9}}=\frac{(J-J^{\ast })}{J^{3}}-\frac{1}{\bar{\unicode[STIX]{x1D6E4}}}\left(\frac{\unicode[STIX]{x1D70E}_{c0}^{2}+2s^{\ast }\unicode[STIX]{x1D70E}_{c0}\unicode[STIX]{x1D70E}_{o0}\cosh (J\,\unicode[STIX]{x1D712}^{\prime })+s^{\ast \,2}\unicode[STIX]{x1D70E}_{o0}^{2}}{\sinh ^{2}(J\,\unicode[STIX]{x1D712}^{\prime })}\right)_{\unicode[STIX]{x1D6F9}},\end{eqnarray}$$

where $J^{\ast }$ is the value of $J$ in the thin film region, and $\unicode[STIX]{x1D712}^{\prime }\equiv (3Ca_{B})^{2/3}/\unicode[STIX]{x1D716}_{2}$ (the ratio between the characteristic or ‘nominal’ film thickness according to Bretherton (Reference Bretherton1961) and the Debye screening length (Wright Reference Wright2007)).

Since the capillary number has been scaled out, we expect $J^{\ast }$ just to be a function of the remaining parameters $\bar{\unicode[STIX]{x1D6E4}}$ , $\unicode[STIX]{x1D70E}_{c0}$ , $\unicode[STIX]{x1D70E}_{o0}$ , $s^{\ast }$ and $\unicode[STIX]{x1D712}^{\prime }$ . The issue however is that $J^{\ast }$ is a priori unknown and must be solved as part of the solution of the problem. This is addressed by rescaling yet again, introducing variables $G=J/J^{\ast }$ and $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6F9}/J^{\ast }$ . We also define

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}^{\ast }(\bar{\unicode[STIX]{x1D6E4}},J^{\ast })=\bar{\unicode[STIX]{x1D6E4}}/J^{\ast } & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D712}^{\ast }(\unicode[STIX]{x1D712}^{\prime },J^{\ast })=\unicode[STIX]{x1D712}^{\prime }J^{\ast }. & \displaystyle\end{eqnarray}$$

The resulting equation is

(2.5) $$\begin{eqnarray}G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}=\frac{(G-1)}{G^{3}}-\frac{1}{\unicode[STIX]{x1D6E4}^{\ast }}\left(\frac{\unicode[STIX]{x1D70E}_{c0}^{2}+2s^{\ast }\unicode[STIX]{x1D70E}_{c0}\unicode[STIX]{x1D70E}_{o0}\cosh (G\,\unicode[STIX]{x1D712}^{\ast })+s^{\ast \,2}\unicode[STIX]{x1D70E}_{o0}^{2}}{\sinh ^{2}(G\,\unicode[STIX]{x1D712}^{\ast })}\right)_{\unicode[STIX]{x1D701}}\end{eqnarray}$$

the final term in brackets as in (2.1) representing the electro-osmotic tension.

We now know that in the thin film region (for small $\unicode[STIX]{x1D701}$ ) $G\rightarrow 1$ , whereas as we approach the capillary static region (for large $\unicode[STIX]{x1D701}$ ) $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ must approach some asymptotic value (this follows because the right-hand side of (2.5) vanishes for large $G$ ). We denote this large $\unicode[STIX]{x1D701}$ asymptotic value by $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ . Based on the way in which we have defined $G$ , $\unicode[STIX]{x1D701}$ , $J$ and $\unicode[STIX]{x1D6F9}$ in terms of $h$ and $Z$ it follows that $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}=J^{\ast }J_{\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6F9}}=J^{\ast }h_{ZZ}$ . Since $h_{ZZ}$ is expected to approach unity in the capillary static region, knowing $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ is sufficient to give $J^{\ast }$ . Equation (2.5) can therefore be solved numerically (for any given $\unicode[STIX]{x1D6E4}^{\ast }$ , $\unicode[STIX]{x1D712}^{\ast }$ , $\unicode[STIX]{x1D70E}_{c0}$ , $\unicode[STIX]{x1D70E}_{o0}$ and $s^{\ast }$ ) and hence the corresponding value of $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ can be determined.

In the Bretherton case, finding $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ completes the solution of the problem. In the presence of electro-osmotic terms, however this is not the full solution to the problem, because $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ has been computed for values $\unicode[STIX]{x1D6E4}^{\ast }$ and $\unicode[STIX]{x1D712}^{\ast }$ but these are not given a priori, but rather are themselves functions of $J^{\ast }$ according to (2.3)–(2.4). We are therefore required to vary a guessed value for $J^{\ast }$ (for fixed $\bar{\unicode[STIX]{x1D6E4}}$ and $\unicode[STIX]{x1D712}^{\prime }$ and therefore variable $\unicode[STIX]{x1D6E4}^{\ast }$ and $\unicode[STIX]{x1D712}^{\ast }$ ) and keep computing $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ until the computed $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701},\infty }$ equals the guessed $J^{\ast }$ .

Finding the required $J^{\ast }$ that satisfies that condition is straightforward if we first solve the Bretherton problem with no electro-osmotic effects to obtain a well-defined value $J_{B}^{\ast }$ , and then recall that a disjoining electro-osmotic tension leads to thicker films (and so must have $J^{\ast }>J_{B}^{\ast }$ ) whereas a conjoining electro-osmotic pressure leads to thinner films (and so has $J^{\ast }<J_{B}^{\ast }$ ). Moreover cases with large $\bar{\unicode[STIX]{x1D6E4}}$ and/or large $\unicode[STIX]{x1D712}^{\prime }$ have weak electro-osmotic effects and hence must be close to the Bretherton case, meaning that we can easily work downwards from large $\bar{\unicode[STIX]{x1D6E4}}$ and/or $\unicode[STIX]{x1D712}^{\prime }$ to smaller values of these parameters.

2.3 Parameters $\bar{\unicode[STIX]{x1D6E4}}$ and $\unicode[STIX]{x1D716}_{2}$

Earlier we gave the physical interpretation of $\bar{\unicode[STIX]{x1D6E4}}$ and $\unicode[STIX]{x1D716}_{2}$ but did not give actual formulae.

Specifically $\bar{\unicode[STIX]{x1D6E4}}$ is defined as

(2.6) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6E4}}=2\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D700}_{w}/Rq^{2}s_{c}^{\ast \,2},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is interfacial tension, $\unicode[STIX]{x1D700}_{w}$ is electrical permittivity of the aqueous phase, $R$ is half-thickness of the capillary, $q$ is the electronic charge, $s_{c}^{\ast }$ is the density of adsorption sites on the wall. Note that Wilmott et al. (Reference Wilmott, Breward and Chapman2018) employed a related parameter $\unicode[STIX]{x1D6E4}$ in lieu of $\bar{\unicode[STIX]{x1D6E4}}$ with a slightly different definition, but in the formulation that we present $\bar{\unicode[STIX]{x1D6E4}}$ is the more convenient parameter. Observe (see e.g. (2.1) or (2.5)) that the larger the value of $\bar{\unicode[STIX]{x1D6E4}}$ the less important the electro-osmotic terms. In the context of waterflooding, an increased $\bar{\unicode[STIX]{x1D6E4}}$ could be achieved by decreasing $R$ (although it would be challenging to extract oil from a rock with extremely narrow pores) or by reducing the surface density of adsorption sites (which presumably could be achieved by changing the chemistry of the clay comprising the capillary wall).

The parameter $\unicode[STIX]{x1D716}_{2}$ is the electrical double layer thickness non-dimensionalised by $R$

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D716}_{2}=(\unicode[STIX]{x1D700}_{w}k_{B}T/(2(C^{+}+4C^{2+})q^{2}R^{2}))^{1/2},\end{eqnarray}$$

where $\unicode[STIX]{x1D700}_{w}$ is the electric permittivity of the aqueous phase, $k_{B}$ is Boltzmann’s constant, $T$ is absolute temperature, $q$ is electronic charge, $C^{+}$ is number density of monovalent ions and $C^{2+}$ is number density of divalent ions. However, the parameter that determines whether or not electro-osmotic effects are screened across the aqueous layer is $\unicode[STIX]{x1D712}^{\prime }$ (itself introduced in (2.2)) which evaluates to

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D712}^{\prime }\equiv (3Ca_{B})^{2/3}/\unicode[STIX]{x1D716}_{2}=(3Ca_{B})^{2/3}(2(C^{+}+4C^{2+})q^{2}R^{2}/\unicode[STIX]{x1D700}_{w}k_{B}T)^{1/2}.\end{eqnarray}$$

Note that Wilmott et al. (Reference Wilmott, Breward and Chapman2018) employed a related parameter $\unicode[STIX]{x1D712}$ in lieu of $\unicode[STIX]{x1D712}^{\prime }$ with a slightly different definition: in the formulation presented here, using the parameter $\unicode[STIX]{x1D712}^{\prime }$ turns out to be more convenient. Larger $\unicode[STIX]{x1D712}^{\prime }$ implies more effective charge screening making electro-osmotic terms less relevant. Practical ways of altering $\unicode[STIX]{x1D712}^{\prime }$ when recovering oil from a given formation include injecting fluid at a different speed (in effect changing the capillary number) or varying the salinity (changing $C^{+}$ and/or  $C^{2+}$ ). In typical systems of interest, the number density of monovalent ions $C^{+}$ tends to be rather larger than that of divalent ions $C^{2+}$ meaning that to a good approximation $\unicode[STIX]{x1D712}^{\prime }\propto (C^{+})^{1/2}$ . As will be discussed below, changing $C^{+}$ and/or $C^{2+}$ not only affects $\unicode[STIX]{x1D712}^{\prime }$ but also has a side effect of changing the charge states of the wall and oil surfaces, denoted $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ . Since the electro-osmotic term in (2.2) is sensitive to all these parameters ( $\unicode[STIX]{x1D712}^{\prime }$ , $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ ), the net effect of changing salinity can potentially be quite complex.

2.4 Determining surface charges on the wall and oil $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$

So far we have stated that the charge per adsorption site on the capillary wall $\unicode[STIX]{x1D70E}_{c0}$ and the oil $\unicode[STIX]{x1D70E}_{o0}$ can vary from $-1$ (no positive ions adsorbed whatsoever) to $+1$ (divalent positive ions adsorbed on all sites). In a number of our calculations we will, for simplicity, impose fixed values for $\unicode[STIX]{x1D70E}_{c0}$ and/or $\unicode[STIX]{x1D70E}_{o0}$ between these extremal values $\pm 1$ . Nevertheless we wish to establish how $\unicode[STIX]{x1D70E}_{c0}$ and/or $\unicode[STIX]{x1D70E}_{o0}$ might be determined more generally.

Following Wilmott et al. (Reference Wilmott, Breward and Chapman2018) we assume that ion adsorption is governed by a Langmuir isotherm. We suppose as above that $C^{+}$ and $C^{2+}$ are monovalent and divalent ion number densities. Moreover we suppose that there are Langmuir constants: $K_{c}^{1}$ (governing monovalent adsorption on the capillary wall), $K_{c}^{2}$ (governing divalent adsorption on the wall), $K_{o}^{1}$ (governing monovalent adsorption on oil) and $K_{o}^{2}$ (governing divalent adsorption on oil). We next define dimensionless ion concentrations $c^{+}$ and $c^{2+}$ by multiplying each ion number density through by $K_{c}^{1}$ , and dimensionless Langmuir constants ( ${\mathcal{K}}_{c}^{1}$ , ${\mathcal{K}}_{c}^{2}$ , ${\mathcal{K}}_{o}^{1}$ and ${\mathcal{K}}_{o}^{2}$ ) by dividing each Langmuir constant through by $K_{c}^{1}$ : by definition ${\mathcal{K}}_{c}^{1}$ equals unity. Note that this differs from the non-dimensionalisation procedure of Wilmott et al. (Reference Wilmott, Breward and Chapman2018) who instead made concentration dimensionless based on an arbitrary concentration that they selected, considered to be typical of a low salinity waterflood.

Dimensionless ion concentrations $c^{+}$ and $c^{2+}$ are of course affected by the electrostatic potential, with concentrations being lower in regions of positive potential, and higher in regions of negative potential. The concentrations of negatively charged balancing counter ions are affected in the opposite way. Throughout this study, as is standard in Debye–Hückel theory (Wright Reference Wright2007), we will assume that the perturbations in ion concentration induced by electrostatic effects are small compared to the unperturbed concentrations $c_{0}^{+}$ and $c_{0}^{2+}$ in the absence of electrostatics. Under this assumption then, these unperturbed concentrations $c_{0}^{+}$ and $c_{0}^{2+}$ will determine the amount of adsorption on the wall and oil surfaces, despite those surfaces having a non-zero electrostatic potential.

In the notation of Wilmott et al. (Reference Wilmott, Breward and Chapman2018), the fraction of adsorption sites on the wall occupied by adsorbed monovalent and divalent positive ions are respectively denoted $s_{c0}$ and $s_{c0}^{+}$ , the analogous fractions for adsorption on oil being $s_{o0}$ and $s_{o0}^{+}$ . In terms of the dimensionless concentrations and dimensionless Langmuir parameters these evaluate to

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle s_{c0}={\mathcal{K}}_{c}^{1}c_{0}^{+}/(1+{\mathcal{K}}_{c}^{1}c_{0}^{+}+{\mathcal{K}}_{c}^{2}c_{0}^{2+}) & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle s_{c0}^{+}={\mathcal{K}}_{c}^{2}c_{0}^{2+}/(1+{\mathcal{K}}_{c}^{1}c_{0}^{+}+{\mathcal{K}}_{c}^{2}c_{0}^{2+}) & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle s_{o0}={\mathcal{K}}_{o}^{1}c_{0}^{+}/(1+{\mathcal{K}}_{o}^{1}c_{0}^{+}+{\mathcal{K}}_{o}^{2}c_{0}^{2+}) & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle s_{o0}^{+}={\mathcal{K}}_{o}^{2}c_{0}^{2+}/(1+{\mathcal{K}}_{o}^{1}c_{0}^{+}+{\mathcal{K}}_{o}^{2}c_{0}^{2+}). & \displaystyle\end{eqnarray}$$

Finally the charges per adsorption site are $\unicode[STIX]{x1D70E}_{c0}=2s_{c0}^{+}+s_{c0}-1$ for the wall and $\unicode[STIX]{x1D70E}_{o0}=2s_{o0}^{+}+s_{c0}-1$ for the oil and hence

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{c0}=({\mathcal{K}}_{c}^{2}c_{0}^{2+}-1)/(1+{\mathcal{K}}_{c}^{1}c_{0}^{+}+{\mathcal{K}}_{c}^{2}c_{0}^{2+}) & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{o0}=({\mathcal{K}}_{o}^{2}c_{0}^{2+}-1)/(1+{\mathcal{K}}_{o}^{1}c_{0}^{+}+{\mathcal{K}}_{o}^{2}c_{0}^{2+}). & \displaystyle\end{eqnarray}$$

The subscripts ‘0’ on all these parameters ( $s_{c0}$ , $s_{c0}^{+}$ , $s_{o0}$ , $s_{o0}^{+}$ , $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ ) reflect the fact that they can be evaluated at ‘unperturbed’ concentrations $c_{0}^{+}$ and $c_{0}^{2+}$ instead of $c^{+}$ and $c^{2+}$ . Moreover in the interests of simplicity, it is assumed that sufficient ions are present in the aqueous film overall such that adsorption onto the capillary wall or oil surfaces does not lead to any significant change in ion concentration in the bulk.

To summarise, for a number of calculations we will perform, fixed values of $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ will be imposed. In cases however when fixed values are not imposed, equations (2.13)–(2.14) will be used instead.

2.5 Pressure drop across the droplet front

So far our discussion has focussed on finding how electro-osmotic effects affect the thin film region, specifically finding the parameter $J^{\ast }$ which determines the thickness of that region. Also of interest however is the pressure drop. At leading order the (dimensionless) pressure drop across the capillary static droplet front is unity, since the (dimensionless) curvature radius at the front is likewise unity. It is of interest however to determine how much the pressure drop at the front is perturbed by droplet motion. The technique for doing this was discussed by Wilmott et al. (Reference Wilmott, Breward and Chapman2018). The finding was that the dimensionless pressure drop $\unicode[STIX]{x0394}p$ could be expressed in the form $1+(3Ca_{B})^{2/3}\unicode[STIX]{x0394}p^{\ast }$ where $\unicode[STIX]{x0394}p^{\ast }$ could be obtained via $\unicode[STIX]{x0394}p^{\ast }\equiv \lim _{\unicode[STIX]{x1D701}\rightarrow \infty }(G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}G-G_{\unicode[STIX]{x1D701}}^{2}/2)$ . It follows that not only must $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ have a well-defined limiting value as $\unicode[STIX]{x1D701}\rightarrow \infty$ at the end of the transition region (thereby determining $J^{\ast }$ ) but also $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}G-G_{\unicode[STIX]{x1D701}}^{2}/2$ must have a well defined limiting value that determines $\unicode[STIX]{x0394}p^{\ast }$ .

This completes our presentation of the model and governing equations. To summarise the key equation we must solve is (2.5) for $G$ versus $\unicode[STIX]{x1D701}$ . The large $\unicode[STIX]{x1D701}$ limiting behaviour of this function $G$ furnishes the information we require, specifically the values of thin film thickness $J^{\ast }$ and the pressure drop change $\unicode[STIX]{x0394}p^{\ast }$ . Parameter values $\unicode[STIX]{x1D6E4}^{\ast }$ , $\unicode[STIX]{x1D712}^{\ast }$ , $s^{\ast }$ , $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ must be supplied (the first two parameters mentioned being related to $\bar{\unicode[STIX]{x1D6E4}}$ and $\unicode[STIX]{x1D712}^{\prime }$ via (2.3)–(2.4), and the final two being related, if required, to Langmuir parameters ( ${\mathcal{K}}_{c}^{1}$ , ${\mathcal{K}}_{c}^{2}$ , ${\mathcal{K}}_{o}^{1}$ , ${\mathcal{K}}_{o}^{2}$ ) via (2.13)–(2.14)): an interpretation of each of these parameters has been summarised in table 1. In the next section the numerical methodology for solving the equations and the choices of parameter values are detailed.

Table 1.

Dimensionless parameters affecting the electro-osmotic term in the model. Note that, although mentioned in the table, ${\mathcal{K}}_{c}^{1}\equiv 1$ by definition here.

3.1 Numerical methodology

We have solved (2.5) using a fourth-order Runge–Kutta method (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1992).

For sufficiently small $\unicode[STIX]{x1D701}$ we know that $G$ should approach unity. Meanwhile for sufficiently large $\unicode[STIX]{x1D701}$ we know that $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ should approach a constant. In order to be able to solve the equation, we need however, in the small $\unicode[STIX]{x1D701}$ limit, an initial condition for both $G_{\unicode[STIX]{x1D701}}$ and $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ , in addition to that for $G$ . Such conditions are obtained by linearising equation (2.5) setting $g\approx G-1$ . It follows

(3.1) $$\begin{eqnarray}g_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}=g-g_{\unicode[STIX]{x1D701}}(\unicode[STIX]{x1D6E4}^{\ast })^{-1}\unicode[STIX]{x2202}T_{EO}/\unicode[STIX]{x2202}G|_{G=1},\end{eqnarray}$$

where $T_{EO}$ is the electro-osmotic tension (i.e. the expression within the brackets in the final term of (2.5)). The solution of (3.1) is $g=\unicode[STIX]{x1D700}\exp (\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D701})$ where $\unicode[STIX]{x1D700}$ is an arbitrarily chosen small parameter (we set $\unicode[STIX]{x1D700}=10^{-5}$ ) and

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}^{3}=1-\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D6E4}^{\ast })^{-1}\unicode[STIX]{x2202}T_{EO}/\unicode[STIX]{x2202}G|_{G=1}.\end{eqnarray}$$

Solving this cubic equation to determine $\unicode[STIX]{x1D6EC}$ is simple. It turns out that there is just one positive real root, and (at the front end of the drop, which is the case we consider here) this is the root we require (Wilmott et al. Reference Wilmott, Breward and Chapman2018). If $\unicode[STIX]{x1D6E4}^{\ast }$ is large, the value of $\unicode[STIX]{x1D6EC}$ is close to unity, whereas if $\unicode[STIX]{x1D6E4}^{\ast }$ is small, the two terms on the right-hand side come into balance (if $\unicode[STIX]{x2202}T_{EO}/\unicode[STIX]{x2202}G|_{G=1}>0$ ) or else the left-hand side comes into balance with the second term on the right (if $\unicode[STIX]{x2202}T_{EO}/\unicode[STIX]{x2202}G|_{G=1}<0$ ). Given an initial estimate for $\unicode[STIX]{x1D6EC}$ we can determine the exact value by Newton–Raphson iteration. At $\unicode[STIX]{x1D701}=0$ we now have $G=1+\unicode[STIX]{x1D700}$ , $G_{\unicode[STIX]{x1D701}}=\unicode[STIX]{x1D700}\,\unicode[STIX]{x1D6EC}$ and $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}=\unicode[STIX]{x1D700}\,\unicode[STIX]{x1D6EC}^{2}$ . We integrated forward via Runge–Kutta with a step size $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}=0.1$ . The choice of step size $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}$ was dictated by the requirement that changes $G$ , $G_{\unicode[STIX]{x1D701}}$ and $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ during any given step should be small compared to the values of $G$ , $G_{\unicode[STIX]{x1D701}}$ and $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ at the start of the step. The integration proceeded until $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ had converged to 6 figures.

An issue however is that the other parameter we seek namely $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}G-G_{\unicode[STIX]{x1D701}}^{2}/2$ converges much more slowly than $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ does. In an attempt to obtain convergence for this parameter we allowed $\unicode[STIX]{x1D701}$ up to very large values (as large as 1000). In that case however both $G_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}G$ and $G_{\unicode[STIX]{x1D701}}^{2}/2$ are then large numbers (of the order of $\unicode[STIX]{x1D701}^{2}$ and hence of the order of a million) and the value we seek is the difference between these two large numbers. To overcome this, we devised techniques to accelerate the convergence, as discussed in appendix C.

3.2 Selecting parameter values: $\unicode[STIX]{x1D712}^{\prime }$ , $\bar{\unicode[STIX]{x1D6E4}}$ and $s^{\ast }$

With the numerical solution technique now specified, it remains to select the sets of parameter values for which solutions will be obtained. Based on the literature, Wilmott et al. (Reference Wilmott, Breward and Chapman2018) proposed a number of parameter values for displacement of an oil droplet from a capillary. Our intention here is to consider parameter values comparable with those of Wilmott et al. (Reference Wilmott, Breward and Chapman2018), but recognising that a number of parameters are sensitive to surface physical chemistry, and these might exhibit quite wide variation depending on the chemical composition of the system (see e.g. Newcombe & Ralston (Reference Newcombe and Ralston1992) for a situation in which surface properties are strongly modified by chemistry).

We estimate a capillary number $Ca_{B}$ of the order of $10^{-7}$ . Obtaining this value requires a conversion between the definition of $Ca$ used by Wilmott et al. (Reference Wilmott, Breward and Chapman2018) (based on an imposed pressure) and the capillary number $Ca_{B}$ used here (based on an imposed injection speed). Since $Ca_{B}$ is proportional to injection speed we estimate that in an oil recovery operation it could vary by an order of magnitude either side of $10^{-7}$ depending on how the operation is realised. Meanwhile, values of the parameter $\unicode[STIX]{x1D716}_{2}$ (the ratio between Debye length and channel half-thickness) cited by Wilmott et al. (Reference Wilmott, Breward and Chapman2018) are between $1.5\times 10^{-3}$ in what they call a low salinity case (corresponding to $17~\text{mol}~\text{m}^{-3}$ of monovalent positive ions) and $2\times 10^{-4}$ (corresponding to $1700~\text{mol}~\text{m}^{-3}$ of monovalent positive ions).

Remembering that $\unicode[STIX]{x1D712}^{\prime }\equiv (3Ca_{B})^{2/3}/\unicode[STIX]{x1D716}_{2}$ , we estimate $\unicode[STIX]{x1D712}^{\prime }\approx 0.03$ in the low salinity case and $\unicode[STIX]{x1D712}^{\prime }\approx 0.22$ in the high salinity case. Given however we commented that $Ca_{B}$ could increase by up to an order of magnitude, and given also that $\unicode[STIX]{x1D712}^{\prime }$ is sensitive to the half-thickness of the pore (see (2.8)) which varies from rock to rock, it is reasonable to suppose that the domain of $\unicode[STIX]{x1D712}^{\prime }$ of interest is somewhere from order unity down to a couple of orders of magnitude smaller than that. To cover this domain specifically we will consider values $\unicode[STIX]{x1D712}^{\prime }$ equal to 5, 2, 1, 0.5, 0.2, 0.1, 0.05 and 0.02. For the $\unicode[STIX]{x1D712}^{\prime }$ values at the upper end of the domain charges are screened and electrostatic effects represent only a small perturbation to the original Bretherton system (Bretherton Reference Bretherton1961). For $\unicode[STIX]{x1D712}^{\prime }$ values at the lower end of the domain however, electrostatic effects should become very significant.

Concerning the parameter $\bar{\unicode[STIX]{x1D6E4}}$ , the estimate is $\bar{\unicode[STIX]{x1D6E4}}\approx 8$ . Obtaining this again requires a conversion from a parameter $\unicode[STIX]{x1D6E4}$ (used by Wilmott et al. (Reference Wilmott, Breward and Chapman2018), defined in terms of a channel length) to the parameter $\bar{\unicode[STIX]{x1D6E4}}$ used here (defined using a channel half-thickness). Since $\bar{\unicode[STIX]{x1D6E4}}$ relies on surface physical chemical properties (in particular it depends on the density of adsorption sites on the capillary wall $s_{c}^{\ast }$ , reported by Li & Xu (Reference Li and Xu2008)) there is likely to be considerable variation in $\bar{\unicode[STIX]{x1D6E4}}$ depending on chemical composition. In view of this, in the interest of simplicity, we set $\bar{\unicode[STIX]{x1D6E4}}=10$ (instead of the original estimated value 8). We also considered a case with an order of magnitude smaller density of adsorption sites and hence considerably weaker electrostatic effects. This case had $\bar{\unicode[STIX]{x1D6E4}}=1000$ . We also considered a case intermediate between these, namely $\bar{\unicode[STIX]{x1D6E4}}=100$ .

Using data from Lewis (Reference Lewis1937), it was estimated (Wilmott et al. Reference Wilmott, Breward and Chapman2018) that there could be 75 times as many adsorption sites on oil as on the wall, i.e.  $s^{\ast }=75$ . Again this estimate is likely to be very sensitive to chemical composition. We chose a value of $s^{\ast }=10$ , significantly smaller than the estimate of Wilmott et al. (Reference Wilmott, Breward and Chapman2018), but still significantly larger than unity (and hence still with more adsorption sites on oil than on the wall). The reason for choosing this value $s^{\ast }=10$ was to explore the effect of interactions between the charged oil and the charged capillary wall. When $s^{\ast }=75$ the charge on the oil completely overwhelms that on the wall. When $s^{\ast }=10$ however the charge on the oil remains dominant, but it becomes possible to detect effects arising from the presence or absence of charge on the wall.

Note one consequence of choosing the combination $\bar{\unicode[STIX]{x1D6E4}}=10$ and $s^{\ast }=10$ . Electro-osmotic effects that refer to the oil (scaling as $s^{\ast \,2}/\bar{\unicode[STIX]{x1D6E4}}$ in (2.2)) and those referring to the oil–wall interaction (scaling as $s^{\ast }/\bar{\unicode[STIX]{x1D6E4}}$ ) can still be important, even though those that refer to the capillary wall (scaling as $1/\bar{\unicode[STIX]{x1D6E4}}$ ) might not be. Furthermore it is possible for the oil–wall interactions to dominate those purely on the oil, as the latter turn out to be more strongly screened at large distances than the former: see appendix B.

3.3 Selecting parameter values: surface charges, Langmuir constants and salinities

Specifying the number of available adsorption sites (via the parameters $\bar{\unicode[STIX]{x1D6E4}}$ and $s^{\ast }$ ) does not completely specify the surface state of either the capillary wall or the oil. It is also necessary to know how many ions actually adsorb. This is determined by the values of Langmuir constants and ion concentrations.

Based on data from Fletcher & Sposito (Reference Fletcher and Sposito1989), it was determined (Wilmott et al. Reference Wilmott, Breward and Chapman2018) that the product of the Langmuir parameter $K_{c}^{1}$ and the low salinity monovalent ion number density $C^{+}$ was approximately $1.7\times 10^{-2}$ . Since a high salinity monovalent ion concentration was taken to be 100 times the low salinity one, and since concentrations are made dimensionless here on the scale $(K_{c}^{1})^{-1}$ , it follows that the dimensionless ion concentrations $c_{0}^{+}$ will vary from roughly order unity (high salinity) down to about two orders of magnitude smaller than that (low salinity).

The work of Wilmott et al. (Reference Wilmott, Breward and Chapman2018) considered a mix of monovalent and divalent ions. In the low salinity system, there were 3.4 times more monovalent ions than divalent ones, whereas in the high salinity system there were 18 times more monovalent ions than divalent, whilst the number of divalent ions also increased 18-fold between the low and high salinity systems. In this work we decided to consider instead 3 different situations in respect of the divalent ions, not identical to the cases of Wilmott et al. (Reference Wilmott, Breward and Chapman2018), but nevertheless in a similar parameter range. We considered a case with a dimensionless concentration $c_{0}^{2+}$ for divalent ions fixed at a low value 0.0025, but with monovalent ion concentration $c_{0}^{+}$ varying between 0.01 and 1.0. We also considered a case with a slightly higher divalent ion concentration $c_{0}^{2+}$ fixed at 0.05, again with monovalent ion concentration $c_{0}^{+}$ varying between 0.01 and 1.0. In addition we considered $c_{0}^{2+}$ : $c_{0}^{+}$ held in a fixed 1:4 ratio, still with $c_{0}^{+}$ varying between 0.01 and 1.0. For each of these cases we also allowed $\unicode[STIX]{x1D712}^{\prime }$ to vary from 0.02 to 0.2 as $c_{0}^{+}$ itself varied between 0.01 and 1 (noting that $\unicode[STIX]{x1D712}^{\prime }$ typically scales like the square root of ion concentration, see § 2.3).

With the salinity levels now specified, we consider the Langmuir parameters. Being surface physical chemistry properties, we anticipate that Langmuir parameters can be very sensitive to chemical composition. Here however we scale the system such that the dimensionless Langmuir constant for monovalent ions on the wall ${\mathcal{K}}_{c}^{1}$ is unity by definition. Based on Fletcher & Sposito (Reference Fletcher and Sposito1989), it was suggested by Wilmott et al. (Reference Wilmott, Breward and Chapman2018) that the analogous constant for divalent ions on the wall ${\mathcal{K}}_{c}^{2}$ could be some 200 times larger than ${\mathcal{K}}_{c}^{1}$ . We select a value with the same order of magnitude but for simplicity we set ${\mathcal{K}}_{c}^{2}=100$ . Regarding the Langmuir constant for monovalent ions adsorbing on oil ${\mathcal{K}}_{o}^{1}$ , Wilmott et al. (Reference Wilmott, Breward and Chapman2018) claimed, based on work of Fournier et al. (Reference Fournier, Oelkers, Gout and Pokrovski1998), that this was just 1.7 times smaller than ${\mathcal{K}}_{c}^{1}$ . To within an order of magnitude, we set for simplicity ${\mathcal{K}}_{o}^{1}={\mathcal{K}}_{c}^{1}$ . Meanwhile ${\mathcal{K}}_{o}^{2}$ was suggested (based on Joseph (Reference Joseph1946)) to be approximately 0.3 times ${\mathcal{K}}_{c}^{1}$ . For simplicity, we chose to set ${\mathcal{K}}_{o}^{2}=0.1$ (an order of magnitude less than ${\mathcal{K}}_{c}^{1}$ or ${\mathcal{K}}_{o}^{1}$ ).

For the above parameter choices, equation (2.14) indicates that the oil is always negatively charged and in fact $\unicode[STIX]{x1D70E}_{o0}$ is well approximated by $-1/(1+{\mathcal{K}}_{o}^{1}c_{0}^{+})$ (which is the result neglecting the effect of the divalent ions). For low salinity systems in particular, with $c_{0}^{+}\ll 1$ , it then turns out that $\unicode[STIX]{x1D70E}_{o0}$ is close to $-1$ , and in fact for simplicity we will perform many of our calculations with $\unicode[STIX]{x1D70E}_{o0}$ set to $-1$ . The situation with $\unicode[STIX]{x1D70E}_{c0}$ is a little more complicated. It is evident from (2.13) that $\unicode[STIX]{x1D70E}_{c0}$ changes sign when $c_{0}^{2+}=1/{\mathcal{K}}_{c}^{2}$ . Since ${\mathcal{K}}_{c}^{2}=100$ here, the sign change occurs when $c^{2+}=0.01$ which is well within the domain of divalent ion concentrations of interest. For simplicity we will consider a number of cases with the charge state of the wall $\unicode[STIX]{x1D70E}_{c0}$ set to either $-1$ , 0 or 1, in addition to cases in which $\unicode[STIX]{x1D70E}_{c0}$ varies according to (2.13).

Changing the sign of $\unicode[STIX]{x1D70E}_{c0}$ may have a relatively limited effect, because (for the current set of Langmuir parameters) the charge on the oil $\unicode[STIX]{x1D70E}_{o0}$ tends to remain negative, and moreover the density of adsorption sites on the oil is rather larger than on the wall (by a factor $s^{\ast }=10$ ). Interesting effects can however be predicted for a different set of Langmuir parameters. If the affinity of the oil for divalent ions is increased significantly such that ${\mathcal{K}}_{o}^{2}=25$ say (with all other Langmuir parameters remaining unchanged), then we could see sign changes for the charges on both the wall and oil, occurring respectively at divalent ion concentrations $c_{0}^{2+}=1/{\mathcal{K}}_{c}^{2}=0.01$ and $c_{0}^{2+}=1/{\mathcal{K}}_{o}^{2}=0.04$ .

For $c_{0}^{2+}$ between these values 0.01 and 0.04, the charges on the wall and oil are opposite in sign. If we assume a $c_{0}^{2+}$ : $c_{0}^{+}$ ratio of 1:20 (comparable with the 1:18 ratio assumed by Wilmott et al. (Reference Wilmott, Breward and Chapman2018) in their high salinity state) then for $c_{0}^{2+}$ very close to 0.035 (and hence $c_{0}^{+}$ very close to 0.7) the charges are opposite and equal, i.e.  $\unicode[STIX]{x1D70E}_{c0}\approx -s^{\ast }\unicode[STIX]{x1D70E}_{o0}$ . In this state, the nature of the electro-osmotic term changes quite dramatically: see appendix B. Rather than having a strong disjoining tension at small distances, we have a conjoining pressure. This dramatic change in behaviour is of potential interest meaning we will analyse this case ${\mathcal{K}}_{o}^{2}=25$ in addition to the original situation with ${\mathcal{K}}_{o}^{2}=0.1$ . This clearly represents a very significant change in the affinity that the oil has for divalent ions, implying a rather different chemical composition compared to our original parameter set, but even so the assumed ${\mathcal{K}}_{o}^{2}=25$ for oil still remains less than ${\mathcal{K}}_{c}^{2}=100$ for the wall.

This completes our discussion of the methodology and parameter values to be used in the study. To summarise, our parametric studies will investigate $\unicode[STIX]{x1D712}^{\prime }$ values varying between 5 and 0.02, $\bar{\unicode[STIX]{x1D6E4}}$ values either 10 or 100 or 1000, $s^{\ast }=10$ , as well as $\unicode[STIX]{x1D70E}_{o0}=-1$ (at least in many of our calculations) along with $\unicode[STIX]{x1D70E}_{c0}$ either $-1$ , 0 or 1. In addition to this we will consider ion concentrations $c_{0}^{+}$ varying between 0.01 and 1, with $c_{0}^{2+}$ being either 0.0025 or 0.05, or else $c_{0}^{2+}:c_{0}^{+}$ being in a $1:4$ or $1:20$ ratio. In such cases $\unicode[STIX]{x1D712}^{\prime }$ , $\unicode[STIX]{x1D70E}_{c0}$ and $\unicode[STIX]{x1D70E}_{o0}$ will be determined based on the ion concentrations, with Langmuir parameters being specified, whilst $\bar{\unicode[STIX]{x1D6E4}}$ and $s^{\ast }$ remain at $\bar{\unicode[STIX]{x1D6E4}}=10$ and $s^{\ast }=10$ . Results of these parametric studies are considered in the next section.

Figure 2.

Thickness of the aqueous film $J^{\ast }$ as a function of $\unicode[STIX]{x1D712}^{\prime }$ (the ratio between the ‘nominal’ thickness of the aqueous film and the Debye length) for various charge states: a singly charged surface ( $\unicode[STIX]{x1D70E}_{c0}=0$ , $\unicode[STIX]{x1D70E}_{o0}=-1$ , $s^{\ast }=10$ ), same sign charges ( $\unicode[STIX]{x1D70E}_{c0}=-1$ , $\unicode[STIX]{x1D70E}_{o0}=-1$ , $s^{\ast }=10$ ) and opposite but unequal charges ( $\unicode[STIX]{x1D70E}_{c0}=1$ , $\unicode[STIX]{x1D70E}_{o0}=-1$ , $s^{\ast }=10$ ). The electro-capillary parameter is $\bar{\unicode[STIX]{x1D6E4}}=10$ . The horizontal line indicates the Bretherton case (without electro-osmotic effects). A line proportional to $1/\unicode[STIX]{x1D712}^{\prime }$ is also shown to guide the eye.