Primary cementing of vertical wells: displacement and dispersion effects in narrow eccentric annuli

Abstract

Laminar miscible displacement flows in narrow, vertical, eccentric annuli show a wide range of interesting physical phenomena, e.g. static layers and channels, dispersive spikes/fronts and various instabilities. These are of relevance to the primary cementing of wells, where gas leakage and greenhouse gas emissions can result from ineffective cementing. The current popular way to model this process is via a Hele-Shaw approach, reducing the Navier–Stokes equations via scaling arguments and then averaging across the annular gap: the two-dimensional (2-D) gap averaged (2DGA) model. This leads to a reduced model that is computationally efficient and represents some important features of cementing flows, but is also deficient in capturing the effects of dispersion. While three-dimensional (3-D) simulations are able to capture dispersion and most other physical phenomena, they suffer from excessive computational times due to the extreme aspect ratios of cementing geometries and non-Newtonian properties of the fluids. Here we extend the 2DGA approach by modelling the high Péclet number limit of miscible duct displacement flows, to take into account the leading-order effects of dispersion. Our modified dispersive 2DGA model is well able to approximate dispersive effects that are observed by gap-averaging the results of 3-D computations, including geometric effects of eccentricity and instabilities. This represents a marked improvement over the 2DGA approach, while keeping the computational advantages of a 2-D model.

1. Introduction

There are many millions of oil and gas wells worldwide, both producing and decommissioned. A key operation in the construction of these wells, as well as deep water wells, geothermal wells, waste storage wells, etc., is primary cementing. Primary cementing is the process of placing cement in the annular space between the steel casing and the newly exposed geological formations. A complete and durable layer of cement that seals hydraulically, is the foremost goal of the cement job. Failure to seal the well allows gas to migrate through the cemented layer to surface, called surface casing vent flow or annular casing pressure, or to migrate away from the well (gas migration), see figure 1. The gas may come from the target reservoir or an intermediate zone. There is acute current interest in such leakage due to greenhouse gas emission implications, both environmentally and politically. Gas leakage is typically ${\rm CH}_4$, but sometimes also ${\rm H}_2{\rm S}$, which has additional health consequences. The ${\rm CO}_2$ storage reservoirs are also accessed via wells that must not to leak. Groundwater contamination and damage to subsurface ecology are two other concerns, also present when the leakage involves oil seepage to the surface. Such occurrences motivate better understanding of all cementing processes.

Figure 1. Well cementing: (a) schematic of fluid sequence and flow path; (b) gas invasion, hazards and common leakage pathways for surface casing vent flow (SCVF) or annular casing pressure (ACP) and gas migration (GM).

Trudel et al. (2019) find that 28.5 % of wells, drilled from 2010 to 2018 in British Columbia, Canada, report SCVF. Baillie et al. (2019) report 28 %–32 % of sites in Southeastern Saskatchewan, Canada, emitting gas. Estimates of mean leakage rates per day vary from $0.5\unicode{x2013}59\ {\rm m}^3$ ${\rm CH}_4\,{\rm day}^{-1}$, per site. These include well-site measurements (bubble tests) and truck-based mobile surveys. Undoubtedly, similar incidence and leakage rates exist in other parts of the world, but regulations are different as is access to well data.

Defects in well integrity are generally either caused during primary cementing or occur in the months/years afterwards. The latter include shrinkage and other chemical effects of prolonged hydration reactions, geomechanically induced cracking (or worse), brine penetration and eventually casing corrosion. The former causes are generally related to the fluid mechanics of the cementing operation, in which the in situ drilling mud is displaced by pumping a sequence of fluids around the well: down inside the casing and back upwards in the annulus. There is very little downhole measurement and evaluation during cementing. Consequently modelling, simulation and lab-scale experiments have become important tools.

Studies have appeared in both technical publications/conferences and in the fluid mechanics literature. In the 1960s–1990s, computational modelling was largely limited to simple hydraulics-style calculations. Thus, one-dimensional (1-D) analyses and the development of hydraulics-based estimates and rule-based design systems evolved (McLean, Manry & Whitaker 1967; Smith 1986; Couturier et al. 1990; Ryan, Kellingray & Lockyear 1992). These recommendations were largely sensible for vertical well cementing, e.g. maintaining density and frictional pressure hierarchies between successive fluids pumped, but were also often conservative. From the early 1990s, model-based 1-D simulators were used with increasing regularity in cementing designs and industry case studies.

A two-dimensional (2-D) model for laminar fluid–fluid displacement flows, based on a Hele-Shaw approach, was introduced first by Bittleston, Ferguson & Frigaard (2002). The key assumption of the model is that the annular gap is much narrower compared with the mean circumference as well as the axial length scale. Thus, the velocity and fluid concentration profiles in the radial direction were radially averaged, which reduces the problem to a 2-D one. This 2-D gap-averaged (2DGA) model was analysed theoretically and numerically by Pelipenko & Frigaard (2004a,b). An iterative augmented Lagrangian algorithm coupled with numerical solution of the concentration equation was applied to produce 2-D simulations. The importance of this transition to the 2DGA framework cannot be overstated. It was now possible to simulate and visualise steady-state annular displacement flows, as well unsteady and unstable flows. The computational methods of Pelipenko & Frigaard (2004b) definitively computed static channels of mud on the narrow side of eccentric annuli, as postulated many years earlier (McLean et al. 1967). The links between the 2DGA and earlier hydraulics-based models were analysed, showing their conservative nature (Pelipenko & Frigaard 2004c).

Using the same framework, the 2DGA model has been used to study horizontal annuli displacement flows (Carrasco-Teja et al. 2008), displacement flows with a moving inner cylinder (Carrasco-Teja & Frigaard 2009, 2010; Tardy & Bittleston 2015)), turbulent and mixed regime displacement flows (Maleki & Frigaard 2017, 2018, 2019). In the past decade these simulations have become an industry standard, often matching/predicting defects measured later in field case studies, and are coded into proprietary software of a number of companies. Although the 2DGA models have proven very successful industrially, two issues have arisen.

Firstly, there have been criticisms due to the simplifications inherent in the derivation. Strictly speaking, the model is a narrow-gap Hele-Shaw model, in which curvature is partly neglected. Some modified versions of the modelling approach include such terms (Tardy & Bittleston 2015). While this will produce slightly different results, it is not an approach that we follow. Curvature effects are only one of many effects neglected in reduced models of this nature. Selecting only some terms of equal order is mathematically inconsistent and not guaranteed to improve overall model accuracy.

A more valid concern is that the approach adopted by Bittleston et al. (2002) over-simplifies the averaging process, where transport of fluid concentrations is concerned. Essentially, it is assumed that the average of the product (of velocity and fluid concentration), is equal to the product of the averages. For turbulent flows this is reasonable. However, for laminar flows, cementing displacements are high Péclet number and far from the Taylor dispersion limit. Instead, advective dispersion dominates and is not accounted for in these models. In Renteria & Frigaard (2020), comparisons were made with displacement experiments in horizontal eccentric annuli. While the 2DGA models were able to predict qualitative behaviours quite well, dispersive and inertial effects were beyond the ability of the 2DGA model to predict. Zhang & Frigaard (2021) have made preliminary comparisons of three-dimensional (3-D) simulations with 2DGA model results. Although the results were consistent with each other regarding aspects such as predicting interface behaviour (steady versus unsteady), there were many interesting phenomena, e.g. static layers, dispersive spikes and instabilities, that were predicted to a lesser degree or not at all.

The second concern is that while 20 years ago it was unfeasible to conduct 3-D simulations of non-Newtonian displacement flows, this is no longer the case. Many interesting and practically relevant phenomena occur on the scale of the annular gap, e.g. wall layers, instabilities, washouts. Others are connected with azimuthal secondary flows, and lastly, we have effects that vary along the annulus over tens of metres, (well inclination, eccentricity, geology). Capturing the small scales while simulating a full wellbore of thousands of metres is still not feasible (in the context of performing repeated calculations for design purposes), but simulating 10–50 m is feasible with modest parallel machines.

A number of recent studies have looked at fully 3-D flows in the cementing context (Kragset & Skadsem 2018; Etrati & Frigaard 2019; Skadsem et al. 2019a; Skadsem, Kragset & Sørbø 2019b), and in particular for studying local irregularity effects. These occur when, for example, the surrounding formation is poorly consolidated and/or the drilling operation has anomalies, leading to washouts, i.e. localised cavities. These features are less effectively studied using 2DGA models, where the length scale assumptions fail. More recently, 2DGA models were used to study systematic large-scale irregularities, for example, slow axial variations in ellipticity and eccentricity (Renteria et al. 2022), while 3-D models were used to capture small scale irregularities (Sarmadi et al. 2022). Returning to the question of dispersion, Sarmadi, Renteria & Frigaard (2021) found very good agreement between 3-D simulation results and the earlier experiments of horizontal annulus displacement flows (Renteria & Frigaard 2020). Simulations also reveal a complex flow structure that is simply not visible in the experiments.

Thus, the concern is how to reconcile 3-D simulations that convincingly represent complex physical phenomena, but are not economical due to the large aspect ratio of the computations, with the need for better engineering designs of cementing flows. The numerical methods also require a high level of knowledge/expertise to ensure good computational performance, are run on parallel clusters, require large meshes and consequently generate large data sets for each run. These are not yet suitable for frequent iterative usage in process design, i.e. run quickly on desktop machines. On the other hand, the 2DGA models are clearly missing phenomena of interest: dispersion and inertial effects. While the neglect of inertia is inherent in Hele-Shaw-style models, we believe that dispersion can be better accounted for. This is the goal of our paper.

An outline is as follows. We firstly introduce an overview of the models used and the scope of the study in § 2. Mass conservation properties are validated and compared for both models in § 2.4. In § 3, we give representative examples of different displacement flow regimes in a vertical eccentric annulus, computed by 2DGA and 3-D models, respectively. These cases show both the successes and failures of the Hele-Shaw approach. Then, we rederive the 2DGA model in § 4, starting from the reduced shear flow equations. A dispersive 2DGA model (D2DGA) is obtained by applying a different assumption regarding the distribution of fluids across the gap. Finally, detailed comparisons between the D2DGA and 3-D model are given in § 5. We find significant improvement of the D2DGA model in capturing the gap-scale dispersion. The paper ends with conclusions and future work (§ 6).

In an eccentric annular displacement flow there are at least 12 dimensionless parameters. The main contribution of the paper is in the derivation of the D2DGA model and demonstration of its effectiveness in capturing gap-scale dispersion. This success rests on the multilayer approach in describing the velocity field variations across the annular gap, and should be equally effective for velocity fields that arise from the usual non-Newtonian models (that account for many of the dimensionless parameters). Thus, our model derivation is kept as general as possible. However, for simplicity in illustrating the results, we only present results for two Newtonian fluids in a vertical annulus in all our examples.

2. Methodology

As discussed in § 1 the main objectives of this work are to expose deficiencies in 2-D gap-averaged displacement models of primary cementing and to seek a viable way to improve them. As always in practical fluid mechanics, the choice of model is a balance between accuracy and computational cost. Here we introduce the different models used in our study.

2.1. The 2DGA model

We outline here the current 2DGA model for annular displacement flows, with more detail later in § 4.1, where we improve the model. We adopt the formulation of Maleki & Frigaard (2017), where the model is more comprehensive than that earlier (Bittleston et al. 2002). The geometry of the narrow eccentric annulus is mapped and unwrapped into a Hele-Shaw cell, see figure 2. The key assumption of the model is that the mean annular gap ($2\hat {d}$) is much narrower compared with both the mean circumference ($2 {\rm \pi}\hat {r}_a^*$), and to the length scale of any axial geometric variations. The narrow gap approximation is that $\hat {d}/({\rm \pi} \hat {r}_a^*) = \delta /{\rm \pi} \ll 1$. The velocity and fluid concentration profiles are then averaged in the radial direction, which reduces the problem to a 2-D one. The length scale assumptions also allow us to assume the flow is locally a shear flow. With the additional assumption that the gap-averaged quantities describe the fluid present, we are able to solve the 1-D shear flow in the direction of the modified pressure gradient locally. This provides the closure expressions that are employed in the 2DGA model. The model consists of a nonlinear elliptic equation for the stream function and a sequence of advection equations for the fluid concentrations.

Figure 2. Geometry of the narrow eccentric annulus mapped to the Hele-Shaw cell, for the 2DGA model.

All dimensional variables are with the ‘hat’ accent, e.g. $\hat {w}_0$, and dimensionless variables without. The governing dimensionless equations are listed as follows in (2.1)–(2.3), we discuss only laminar flow regimes in this paper:

(2.1)\begin{gather} \frac{\partial}{\partial t}[Hr_{a}\bar{c}_{k}]+\frac{\partial}{\partial \phi}[H\bar{v}\bar{c}_{k}]+\frac{\partial}{\partial \xi}[r_a H\bar{w}\bar{c}_{k}]=0,\quad k=1,2,\ldots,K, \end{gather}

(2.2)\begin{gather} \boldsymbol{\nabla}_{a} \boldsymbol{\cdot} ( 2r_a H\bar{v},2r_a H\bar{w}) = 0, \quad \Rightarrow\quad \frac{\partial \varPsi }{\partial \phi }=2Hr_{a}\bar{w}, \quad \frac{\partial \varPsi }{\partial \xi }={-}2H\bar{v}, \end{gather}

(2.3)\begin{gather} {\boldsymbol{\nabla}_{a}\boldsymbol{\cdot} [ \boldsymbol{S} + \boldsymbol{b}] = 0.} \end{gather}

Equation (2.1) is a transport equation for the fluid concentrations $\bar {c}_{k}$ for $k=1,2,\ldots,K$, which are advected with the gap-averaged velocity $(\bar {v},\bar {w})$. The scaled half-gap-width is $H(\phi,\xi )$ and the mean radius at depth $\xi$ is $r_a(\xi )$. In laminar miscible flows the effects of molecular diffusion are minimal on the time scales and length scales of relevance. Thus, diffusive effects are absent on the right-hand side of (2.1), but not for turbulent regimes, where Taylor dispersion dominates (Maleki & Frigaard 2017). In deriving (2.1), azimuthal and axial lengths have been scaled with a representative (half-)circumference of the annulus ${\rm \pi} \hat {r}_a^*$. Velocities have been scaled with a mean pumping velocity $\hat {w}_0$, and times with ${\rm \pi} \hat {r}_a^*/\hat {w}_0$. The mean radius $\hat {r}_a^*$ is defined by averaging along the annular flow path. The gap averaged velocity field is incompressible, allowing it to be expressed in terms of a stream function $\varPsi$, i.e. (2.2). The operator ${\nabla }_{a}$ is a radial divergence operator, i.e.

(2.4)\begin{equation} {\nabla}_{a} = \left( \frac{1}{r_{a}}\frac{\partial }{\partial \phi }, \frac{\partial }{\partial \xi } \right), \end{equation}

with ${\nabla }_{a}$ is the corresponding gradient operator.

The stream function satisfies the elliptic equation (2.3), where

(2.5a,b)\begin{equation} \boldsymbol{S} = \left[\frac{r_{a}\tau_w (\left | {\nabla}_{a} \varPsi \right |)}{H \left | {\nabla}_{a}\varPsi \right |} \right] {\nabla}_{a}\varPsi ,\quad \boldsymbol{b} = \frac{r_a (\rho-1)}{{F^2}} (\cos \beta, \sin {\rm \pi}\phi \sin \beta ) . \end{equation}

The function $\boldsymbol {S}$ represents a modified pressure gradient, which we see has magnitude $r_a \tau _w ( | {\nabla }_{a} \varPsi |)/H$. Here $\tau _w$ is the wall shear stress. These quantities have been made dimensionless using a stress scale $\hat {\tau }_0$. In the definition of $\boldsymbol {b}$, $\rho$ is the density (scaled with the density of the first fluid), and $F$ is the Froude number,

(2.6)\begin{equation} {F} = \sqrt{\frac{\hat{\tau}_{0}}{ \hat{\rho}_1 \hat{g} \delta_0 \hat{r}_{a}^* }} . \end{equation}

For clarity, here we fix $\hat {\tau }_{0}=\hat {\mu }_1 \hat {w}_0/\hat {d}$, i.e. based on a viscous stress in the displaced fluid. However, more generally $\hat {\tau }_{0}$ could be chosen to represent any relevant stress (Maleki & Frigaard 2017), depending on the dominant flow regime and rheology. Thus, we see that $\boldsymbol {b}$ represents the effects of buoyancy in this model; together with the imposed flow, buoyancy gradients strongly influence the flow.

The dependency of $\tau _w$ on the local flow rate (mean velocity), on the fluid concentrations and on the rheological parameters will depend on the specific fluid properties, i.e. the model is completed with a closure expression for $\tau _w$. Generally speaking, in laminar flow regimes cementing involves shear thinning, yield-stress fluids. The yield stress of these fluids $\tau _Y$ plays an important role industrially, defining a limiting pressure gradient $r_{a}\tau _{Y}/H$, that must be exceeded locally in order for the fluids to flow; otherwise $| {\nabla }_{a} \varPsi |=0$ and the fluids are stuck. The 2-D model was first analysed by Pelipenko & Frigaard (2004a,b) and there are a large number of works that use this underlying framework.

Numerically, we follow the procedures described in Maleki & Frigaard (2017). The equations are discretised on a regular rectangular mesh, with staggered mesh for the stream function versus concentration. We see that the flow evolves in time only via (2.1), which is time-advanced explicitly using a flux corrected transport scheme. At each time step the nonlinear elliptic stream function equation (2.3), is solved using an iterative augmented Lagrangian algorithm. The velocities are constructed from $\varPsi$ via (2.2) and used for the fluid concentrations in (2.1).

2.2. The 3-D model

The 3-D numerical simulations are computed using a finite volume method within the open source computational fluid dynamics toolbox OpenFOAM (http://www.openfoam.com). The code used here is an updated version of that used in Etrati & Frigaard (2018) for studying miscible displacement flows in pipe flows, wherein the results were successfully compared with experiments. More recently, the code has been used by Sarmadi et al. (2021), comparing with the experimental results of Renteria & Frigaard (2020). Briefly, the full Navier–Stokes equations are solved, coupled to an advection equation for the concentration of fluid 2, which is the displacing fluid. The volume of fluid method is used to capture the interface between the fluids. The momentum, continuity and concentration equations are

(2.7)\begin{gather} \hat{\rho}\left(\frac{\partial \hat{\boldsymbol{u}}}{\partial \hat{t}}+\hat{\boldsymbol{u}}\boldsymbol{\cdot} \hat{\nabla} \hat{\boldsymbol{u}} \right)={-}\hat{\nabla} \hat{p}+\hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\hat{\boldsymbol{\tau}}+\hat{\rho}\hat{\boldsymbol{g}}, \end{gather}

(2.8)\begin{gather} \hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\hat{\boldsymbol{u}}=0, \end{gather}

(2.9)\begin{gather} \frac{\partial c}{\partial \hat{t}}+\hat{\boldsymbol{u}} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} c=0. \end{gather}

Here, $\hat {\boldsymbol {u}}$ and $\hat {p}$ are the velocity and pressure; $\hat {\boldsymbol {\tau } }=\hat {\mu }\hat {\dot {\boldsymbol {\gamma } }}(\hat {\boldsymbol {u}})$ is the deviatoric stress tensor, with $\hat {\dot {\boldsymbol {\gamma } }}$ the rate of strain tensor. In this paper we only consider two Newtonian fluids. The concentration $c\in [0,1]$ enters via the fluid properties: $\hat {\rho }$ and $\hat {\mu }$ are the mixture density and viscosity, defined by

(2.10)\begin{gather} \hat{\rho} =(1-c)\hat{\rho}_{1}+c\hat{\rho}_{2}, \end{gather}

(2.11)\begin{gather} \hat{\mu}=\hat{\mu} _{1}\exp\left[c\,\mathrm{log}\left(\frac{\hat{\mu} _{2}}{\hat{\mu} _{1}}\right)\right]. \end{gather}

The subscript ‘1’ represents the displaced fluid and ‘2’ represents the displacing fluid. We have also used linear interpolation for the viscosity closure, with only mild differences.

A uniform velocity is imposed at the inlet, no-slip is applied at the annulus walls and an outflow condition set at the outlet. The pressure is fixed at the outlet. The PIMPLE algorithm is used to solve the momentum equation, and the implicit second-order Crank–Nicolson method is employed for the time discretisation. Mesh sensitivity/convergence results are given in Sarmadi et al. (2021). Figure 3 illustrates the computational mesh used for one of our cases with eccentricity $e = 0.4$, with $20 \times 100 \times 400$ cells in the radial, azimuthal and axial directions. Figure 3(a) shows half of the annulus, taken from the midplane ($x = 0$), and figure 3(b) shows a middle slice surface. The mesh is refined close to the wall to capture residual wall layers, as common in displacement flows. The computations are carried out in parallel on a multiprocessor machine, generally using 24 cores.

Figure 3. Mesh resolution in (a) axial (vertical) direction and (b) the cross-section for an eccentric case with $e = 0.4$.

2.3. Scope of the study

The aim is to study and improve the 2DGA model, making use of the 3-D computations. We also later wish to use these same models for the purposes of experimental comparison. Thus, we set the annulus dimensions and other parameters of our study based on the experimental set-up of Renteria & Frigaard (2020), which we have also used for our experimental study. The total length of the annulus is 4.8 m, the outer and inner radii are 22.23 mm and 17.46 mm, respectively. This gives a scaled length $Z \approx 62$ for the 2DGA model, i.e. the length is $\approx 62 \times$ the mean circumferential length scale ${\rm \pi} \hat {r}_a^*$, (and is $>1000 \times$ the mean annular gap width). The eccentricity can be adjusted from 0 to 1. Unlike Renteria & Frigaard (2020) and Sarmadi et al. (2021), here we study vertical annuli.

For two Newtonian fluids the flow depends on at least seven dimensionless parameters, and more if the fluids are non-Newtonian. Having adopted $\hat {\tau }_{0}=\hat {\mu }_1 \hat {w}_0/\hat {d}$, as our (viscous) stress scale, the key dimensionless parameters can be taken as the Reynolds number, $Re=\hat {\rho }_{1}\hat {w}_{0}\hat {d}/\hat {\mu }_{1}$, the viscosity ratio, $m=\hat {\mu } _{1}/\hat {\mu } _{2}$, the buoyancy number, $b=(\hat {\rho }_{2}-\hat {\rho }_{1} )\hat {g}\hat {d}^{2}/(\hat {\mu }_{1}\hat {w}_{0})$, a density ratio (or Atwood number), the pipe inclination and eccentricity. We see that $b$ represents the size of the buoyancy vector $\boldsymbol {b}$ in the 2DGA model. Rather than explore all parameters systematically, we instead study a selection of 10 cases that expose key features of vertical displacement. The cases and corresponding parameters are listed in tables 1 and 2. For simplicity, all fluid pairs are Newtonian and the annulus is vertical. Each simulation case is run for a time equivalent to two annular volumes of pumped fluid, using both styles of model.

Table 1. Fluid properties for the designed simulation cases.

Table 2. Corresponding dimensionless parameters for the designed simulation cases.

2.4. Mass conservation properties

To verify that both models are correctly programmed and solved, we select cases 2, 4, 6 and 8, each of which has different flow rates. We integrate the equations forward and compute the total volume of the displacing fluid in the annulus $(\hat {V} _{2})$ at successive time steps. This is plotted in figure 4, in which we compare the two computations with the theoretical volume (from the inflow rate). Stopping the simulations at any particular time and dividing by the annulus volume, gives a fraction of the annulus displaced, i.e. the displacement efficiency. We see that in each case, the 2DGA and 3-D models compute near-identical volumes up until the breakthrough time, i.e. that time at which the displacing fluid first exits the annulus. The breakthrough time is generally earlier for the 3-D simulation, due to dispersion ahead of the main front. We also observe earlier breakthrough for cases 6 and 8, for both of which the displacing fluid is less viscous than the displaced fluid ($m > 1$). The same verification is made for the other cases, with similar agreement; not shown for reasons of brevity.

Figure 4. Verification of volume conservation for the 2DGA and 3-D models: (a) cases 2 and 4; (b) cases 6 and 8. The black lines represent the total annulus volume. The solid lines with other colours represent the theoretical value based on the constant flow rate. The 2DGA and 3-D results are marked by different symbols, as indicated.

3. Successes and failures of the Hele-Shaw approach

Now we explore some of these cases in more detail, here and in the Appendix. We start with case 4, which is representative of a successful displacement flow, (commonly found if $e$ is modest and $b$ is large). Figure 5 shows the comparison of 2-D gap-averaged concentration of case 4 at successive times as indicated. The variable $\phi \in [0,1]$ measures the azimuthal angle from the wide side (indicated by W in the figure), to the narrow side (indicated by N in the figure). The variable $\hat {\xi }$ measures distance along the axis of the well. The thin dashed white line represents the interface between the two fluids, (with averaged concentration $\bar {c}=0.5$). Whereas the 2DGA model computes only the gap-averaged quantities only, the 3-D model computes variables across the annular gap. In order to compare with 2DGA results we post-process by averaging the solution variables from the 3-D model across the annular gap in the radial direction.

Figure 5. Gap-averaged concentration colour maps of case 4 (parameters: $e=0.6$; $Re=20$; $m=0.2$; $b=1000$) computed by (a) 2DGA model; (b) 3-D model.

From the comparison, both computations show a flat interface that evolves steadily along the annulus. Although there is a significant eccentricity the large buoyancy number leads to a steady displacement. In terms of predicting the overall behaviour of the (gap-averaged) 3-D code, the 2DGA is clearly effective. The main difference is in the intensity of the (red coloured) concentration close to the interface. The concentration is significantly less dispersed for the 2DGA model. Accordingly, the front of the displacing fluid of the 3-D model is slightly more advanced than that in figure 5(a). To the eye, it may appear that the right-hand sides of the snapshots of figure 5(b) have concentration that is significantly lower than figure 5(a), making the previous volumetric comparison in figure 4(a) appear doubtful. However, note that the differences in the colour maps are largest on the narrow side of the annulus, where the volume is smallest. Equally, there is a small amount of dispersion ahead of the $\bar {c}=0.5$ contour which is not easily seen. To quantify this dispersive aspect in an averaged way we plot a spatiotemporal map in figure 6. Concentrations from both models are averaged over the cross-section at fixed heights $\hat {\xi }$, to give $\bar {C}(\hat {\xi },\hat {t})$. The linear advance of the front is evident and very slightly faster in the 3-D case. Most of the dispersion is found just behind the front in the 3-D case.

Figure 6. Spatiotemporal colour maps of cross-section-averaged concentration $\bar {C}(\hat {\xi },\hat {t})$ of case 4 computed by (a) 2DGA model; (b) 3-D model.

If figure 5 were typical of all displacement flows, there would be little motivation for improving the 2DGA model. However, we shall see in other examples that although the 2DGA model captures the qualitative features of the 3-D flow, the quantitative differences can be significant. To understand the origin of these differences, we must first ask where the dispersion comes from, and secondly better understand the 2DGA model limitations. Considering the underlying assumptions of the 2DGA model, there are three principal simplifications: (i) neglect of inertial effects in the momentum balance; (ii) elimination of gradients of the viscous stress tensor, simplifying to a shear flow; (iii) simplification of the advective mass transport terms in (2.1). It is the last of these that has a clear dispersive effect.

Regarding the actual dispersion, note that in both models we have no diffusive terms in the concentration equations. For the ranges of typical parameters in these displacement flows, the Péclet numbers would typically be in the range $10^7\unicode{x2013}10^9$, which motivates the neglect of diffusion terms, although they can also be included (Maleki & Frigaard 2017). The numerical methods used here do in fact have well-controlled numerical diffusion, limiting significant smearing of the interface to a few mesh cells. However, intermediate concentrations that are computed at the interface are advected away strongly by secondary flows, regardless of diffusion. In order to effectively represent very small physical diffusivities would require a very fine mesh close to the interface, which is not practically feasible for general studies. Even with minimal diffusion, we would observe dispersion in gap-averaging the 3-D velocity profiles, which are Poiseuille-like.

Although the flow of figure 5(a) has a very sharp interface, illustrating the efficacy of the 2-D numerical scheme in limiting numerical diffusion, there is in fact considerable dispersion even here. To see this, note that with $e=0.6$, the far-field flow is significantly faster on the wide side of the annulus than the narrow side, but yet the displacement front advances steadily. As a rough guide, the ratio of wide to narrow side mean velocities in the far-field is $\approx (1+e)^2/(1-e)^2$ for Newtonian fluids. The steadily propagating front therefore requires a significant secondary azimuthal flow, moving fluid from wide to narrow (and vice versa). This secondary flow is localised close to the advancing front. Thus, considering figure 5(a,b), both models have secondary flows and consequent dispersion, but the effects are different.

Now we consider case 7, which has smaller eccentricity and buoyancy number. Figure 7 shows the gap-averaged concentration colour maps computed by both models. As the light blue colour ahead of the interface is not very clear to see, we plot yellow contour lines which represent $\bar {c}=0.01$. There are distinct differences between models. The 2DGA displacement is steady and piston-like at the front. For the 3-D model the interface ($\bar {c}=0.5$) is near-horizontal on the wide side, but begins to lag behind on the narrow side. Obvious dispersion is seen in the 3-D result. This is again primarily behind the displacement front, but we do see a region extending ahead of the front on the wide side. The dispersion is amplified by the small buoyancy number. Note that here $e=0.2$, but the effects of dispersion along the narrow side are far more significant than in case 4. There is no residual layer in the narrow side of the annulus, as the displacing fluids keep moving. There are also mild instabilities that start close to the interface in the narrow side, as the colour map shows at 85 s and 120 s, which do not apparently originate in the 2DGA physics, even though they are reminiscent of fingering-type instabilities.

Figure 7. Gap-averaged concentration colour maps of case 7 (parameters: $e=0.2$; $Re=20$; $m=0.5$; $b=10$) computed by (a) 2DGA model; (b) 3-D model.

To have a better understanding of the 3-D flow, figure 8 shows snapshots of the isosurface $c=0.2$, at different times for case 7, exploring the annulus section from 2.4 m to 3.6 m, i.e. concentrations of 20 % and above are enclosed within the isosurfaces shown. The colour map on the equally spaced cross-sections represents the dimensional axial velocity. We can clearly see the dispersion, along the centreline of each annular gap and faster on wide side. The dispersive velocity reached two times the imposed mean velocity. The gap-averaged velocity of the narrow side is smaller than that of the wide side, which results in the interface shown in figure 7(b), beginning to extend downstream. In addition, we show an annulus cross-section at $\hat {\xi }=3$ m (the middle of the annular sections of figure 8), to illustrate the azimuthal flow, see figure 9. The colour map shows the concentration distribution and black lines represent streamlines. We see a symmetric flow from narrow side to wide side before the main current passes (figure 9a,b), and the azimuthal flow is not that severe. Then, as the front displaces the secondary flow becomes stronger (figure 9c,d), and appears to lose its symmetry. The flow is mostly from wide side to narrow side, which appears to cause the instabilities observed in figure 7(b). As we gap-averaged the velocity field, the 2DGA model cannot represent the complex asymmetric patterns seen in figure 9(c,d), at least in its current form. The question is, whether such instabilities and asymmetries are caused by inertial or by gap-scale dispersive effects, which are not in the 2DGA model.

Figure 8. Snapshots of the 3-D displacement flow at different times for case 7, (parameters: $e=0.2$; $Re=20$; $m=0.5$; $b=10$): (a) $\hat {t}=55\,{\rm s}$; (b) $\hat {t}=65\,{\rm s}$; (c) $\hat {t}=75\,{\rm s}$; (d) $\hat {t}=85\,{\rm s}$. Colour map of equally spaced slices represents dimensional axial velocity. The wide side of the annulus is on the top.

Figure 9. Azimuthal flow of the interface located at $\hat {\xi }=3$ m at different times for case 7: (a) $\hat {t}=55\,{\rm s}$; (b) $\hat {t}=65\,{\rm s}$; (c) $\hat {t}=75\,{\rm s}$; (d) $\hat {t}=85\,{\rm s}$. The colour map shows the concentration distribution and yellow lines represent streamlines. The wide side of the annulus is on the top.

Further examples of the 2DGA model output and its comparison with 3-D simulations are given in the Appendix. As well as dispersion, these also show significant instabilities (more than those of figure 7), that are not captured by the 2DGA model. Together with the results above, these provide motivation for improving the 2DGA model, so as to better represent cementing displacements. The obvious missing ingredient is to account for dispersion due to the velocity profiles across the annular gap.

4. Accounting for gap-scale dispersion

To account for the gap-scale dispersion in the 2DGA we have two main options. First, we might resolve the flow on the scale of the gap using a simplified 3-D approach. For example, we might follow Tardy et al. (2017) and Tardy (2018) and discretise in the radial direction, while keeping the scaling arguments that lead to a shear flow at leading order. Evaluating the radial variation then allows us to advance the concentration using a velocity field that is 3-D and hence resolve the products of $c$ and the velocity directly. This leads to a model that is somewhere between 2-D and 3-D (say 2.5-D). The additional meshing in the radial direction means larger data sets to handle and slower computation, while still not resolving the inertial effects (as is done in a fully 3-D computation). Nevertheless, with care directed at the computational aspects, the 2.5-D approach of Tardy demonstrably works well for many cases.

The second option is to retain the 2-D description and to derive expressions that can express the gap average of the product of concentration and velocity, in a way that approximates the actual dispersive behaviour. The 2DGA model is compatible with the idea that dispersion is primarily in the direction of the gap-averaged flow. In Maleki & Frigaard (2017) this idea is made more explicit in modelling the Taylor dispersion in turbulent annular flows. However, here the laminar flow is far from the Taylor dispersion regime. Instead, the flows are characterised as large Péclet number, small aspect ratio miscible displacement flows on the annulus gap-scale.

These regimes were studied by Yang & Yortsos (1997), who considered asymptotic solutions to the problem of miscible Newtonian displacement in long channels. Their leading-order transverse flow equilibrium (TFE) model combines a lubrication approximation for the velocity, with a convection–diffusion equation for the fluid concentration, containing the transverse diffusion only. Using this approach, a first-order conservation law was derived for the average concentration of displacing fluid ($\bar {c}$), along the plane channel,

(4.1)\begin{equation} \frac{\partial \bar{c}}{\partial t} + \frac{\partial }{\partial x}F_p = 0, \end{equation}

where $F_p$ is the flux of displacing fluid along the channel. In the limit of zero diffusion $F_p$ can be calculated explicitly and expressed as $F_p(\bar {c};m)$. Solutions of (4.1) exhibit a range of behaviours, including shocks and dispersive spreading. Similar flow regimes were analysed later by Lajeunesse et al. (1999) who extended the theory of Yang & Yortsos (1997) to include buoyancy effects. The propagation is governed again by (4.1), but now the flux is given by

(4.2)\begin{equation} F_p(\bar{c},m,U) = \frac{\bar{c}}{2}\frac{[(2m-3)\bar{c}^2+3]}{[1+(m-1)\bar{c}^3]}+\frac{\bar{c}^2(1-\bar{c})^3}{4U}\frac{[(4m-3)\bar{c}+3]}{[1+(m-1)\bar{c}^3]}. \end{equation}

The first term is identical to that of Yang & Yortsos (1997). Now $F_p$ depends also on additional parameter $U$, as follows:

(4.3)\begin{equation} U = \frac{3{\hat{\mu}_{1}\hat{w}_{0}}}{\hat{h}^{2}\hat{g}(\hat{\rho}_{1}-\hat{\rho}_{2})} .\end{equation}

In the model of Lajeunesse et al. (1999), $\hat {h}$ represents the channel half-width, and fluid 1 is displaced by fluid 2. The flow is vertical, in the direction of gravity and $U \ge 0$, i.e. density stable. In Lajeunesse et al. (1999) there is typographical error for (4.2): corrected in the PhD thesis (Lajeunesse 1999), and above. The expression (4.2) has been analysed in Lajeunesse (1999), to predict flow regimes, which mostly agree with those observed experimentally. The question we now address is whether a similar methodology can be effectively used in our 2-D flows, where the channel width and buoyancy direction varies, and where the fluids may also be non-Newtonian.

4.1. The Hele-Shaw setting

If one considers primary cementing in a vertical casing to represent a Hele-Shaw cell of varying gap, it is tempting to resolve the displacement flow in vertical slices along the length of the well. On each such slice we might expect the above analysis to be valid and deliver a good approximation of the gap-scale dynamics. While correct, this neglects the effects of azimuthal secondary flows which can be significant, as we have seen; fortunately, there is a better solution.

In the derivation of the 2DGA model, (2.1)–(2.3), the Navier–Stokes equations are made dimensionless using a different length scale in the radial direction compared with that in azimuthal and axial directions. This leads to an asymptotic approximation (in terms of the aspect ratio), which to leading order is a shear flow, i.e.

(4.4)$$\begin{gather} 0 ={-} \frac{\partial p}{\partial y}, \end{gather}$$

(4.5)$$\begin{gather}0 ={-} \frac{1}{r_a}\frac{\partial p}{\partial \phi} +\frac{\partial}{\partial y} \tau_{\phi y} + \frac{\rho \sin \beta \sin {\rm \pi}\phi}{{F^2}} , \end{gather}$$

(4.6)$$\begin{gather}0 ={-} \frac{\partial p}{\partial \xi } + \frac{\partial}{\partial y} \tau_{\xi y} - \frac{ \rho \cos \beta }{{F^2}}, \end{gather}$$

where $y \in [-H(\phi,\xi ),H(\phi,\xi )]$ is the scaled radial (gap-)coordinate. These equations have scaled the pressure gradients in order to balance with the leading-order shear flow. The density above has been scaled with the density of fluid 1. We deviate slightly from Maleki & Frigaard (2017) above, since we have chosen to not subtract the static pressure field of fluid 1 from the pressure before scaling. This is because evaluating buoyancy effects requires more care. The leading-order strain rate is

(4.7)\begin{equation} \dot{\gamma} \sim \left[ \left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial w}{\partial y}\right)^2\right]^{1/2}. \end{equation}

Thus, to leading order the constitutive laws depend only on the $y$-derivatives of both $v$ and $w$, and on the fluid present. By integrating across the gap, it is found that the motion is driven in the direction of the modified pressure gradient,

(4.8)\begin{equation} \boldsymbol{G} = (G_{\phi},G_{\xi}) \equiv \left(- \frac{1}{r_a}\frac{\partial p}{\partial \phi} + \frac{ \rho \sin \beta \sin {\rm \pi}\phi}{{F^2}}, -\frac{\partial p}{\partial \xi } - \frac{ \rho \cos \beta }{{F^2}} \right). \end{equation}

In other words, the vector of averaged velocities $(\bar {v},\bar {w})$ is parallel in the $(\phi,\xi )$ plane to $\boldsymbol {G}$.

The difficulty comes in resolving the averaged velocities. In Bittleston et al. (2002) and Maleki & Frigaard (2017), the simplification was made that the concentrations were constant across the annular gap. With this assumption, the density $\rho$ and any rheological properties should be defined as functions of the (mean) fluid concentrations. For the density a simple volume average is appropriate. Having specified the rheology of the mixture, the relevant closure expression can also be derived, based on a plane Poiseuille flow. However, it is clear that these simplifications negate dispersion that will occur due to layering of the fluids across the gap during displacement. Equally, buoyancy effects are limited to gradients of the mean density and not individual fluids.

4.2. Dispersive Hele-Shaw setting

In principle, the same method as Bittleston et al. (2002) and Maleki & Frigaard (2017) can be employed here, but with different assumptions made regarding the distribution of fluids across the gap. For simplicity we assume only two fluids: fluid 1 displaced by fluid 2. We assume that the displacing fluid can disperse (symmetrically) down the centre of the local annular gap. Given the large $Pe$ we assume that fluid 2 occupies the centre of the channel, between interfaces at $y = \mp y_i$; fluid 1 occupies the wall layers: $[-H,-y_i)$ and $(y_i,H]$ (see figure 10). The density $\rho$ in (4.5) and (4.6) will now be constant in each fluid layer, again allowing integration, but now with a jump in fluid properties at $y = \mp y_i$.

Figure 10. Schematic of dispersive Hele-Shaw model setting: (a) narrow annulus, with slice along the local flow direction $\boldsymbol {s}$; (b) assumed configuration of fluid 1 (blue) and fluid 2 (red); local concentration $c = y_i/H$.

The concentration of displacing fluid 2 in the gap we shall denote $\bar {c}(\phi,\xi,t)$, which is given geometrically by $\bar {c} = y_i/H$. The evolution of $\bar {c}$ (neglecting diffusive terms), is governed by

(4.9)\begin{equation} \frac{\partial}{\partial t}[Hr_{a}\bar{c}]+ \boldsymbol{\nabla}_a \boldsymbol{\cdot} \boldsymbol{q} = 0, \end{equation}

where the flux functions are

(4.10)\begin{equation} \boldsymbol{q} = (q_\phi , q_\xi ) = r_a\left( \int^{y_i}_0 v\,{{\rm d} y}, \int^{y_i}_0 w\,{{\rm d} y} \right). \end{equation}

Note that if $(v,w)$ are approximated by their gap-averaged quantities, then (4.9) is exactly (2.1). The questions therefore are, how to evaluate the above fluxes and are we able to specify them in terms of gap-averaged quantities?

The shear stresses in (4.5) and (4.6) are coupled only through the effective viscosity of each fluid, $\eta _1$ and $\eta _2$, which depend on $\dot {\gamma }$. On assuming symmetry about $y=0$, integrating (4.5) and (4.6) for $y \in [0,H]$, and imposing continuity of the shear stresses at $y=y_i$, we find

(4.11)\begin{gather} \frac{\partial v}{\partial y} = \left\{ \begin{array}{@{}ll} \displaystyle{\dfrac{y}{\eta_2(y)}\left( \dfrac{1}{r_a}\dfrac{\partial p}{\partial \phi} - \dfrac{\rho \,f_{\xi}}{r_a} \right) + \dfrac{y}{\eta_2(y)} \dfrac{(H-y_i) \Delta \rho \,f_{\xi} }{H r_a} } , & y \in [0,y_i],\\ \displaystyle{ \dfrac{y}{\eta_1(y)} \left( \dfrac{1}{r_a}\dfrac{\partial p}{\partial \phi} - \dfrac{\rho \,f_{\xi}}{r_a} \right) + \dfrac{(H-y)}{\eta_1(y)} \dfrac{y_i\Delta \rho \,f_{\xi} }{Hr_a} }, & y \in (y_i,H], \end{array} \right. \end{gather}

(4.12)\begin{gather} \frac{\partial w}{\partial y} = \left\{ \begin{array}{@{}ll} \displaystyle{\dfrac{y}{\eta_2(y)}\left(\dfrac{\partial p}{\partial \xi} + \dfrac{\rho \,f_{\phi}}{r_a} \right) - \dfrac{y}{\eta_2(y)} \dfrac{(H-y_i) \Delta \rho \,f_{\phi} }{H r_a} } , & y \in [0,y_i],\\ \displaystyle{ \dfrac{y}{\eta_1(y)} \left(\dfrac{\partial p}{\partial \xi} + \dfrac{\rho \,f_{\phi}}{r_a} \right) - \dfrac{(H-y)}{\eta_1(y)} \dfrac{y_i\Delta \rho \,f_{\phi} }{H r_a} }, & y \in (y_i,H], \end{array} \right. \end{gather}

where $\Delta \rho = \rho _1-\rho _2$ is the density difference, $\rho = (1-y_i/H) \rho _1+ (y_i/H)\rho _2$ is the mean density and

(4.13)\begin{equation} \boldsymbol{f} = (\,f_{\phi},\,f_{\xi}) = \left( \frac{r_a \cos \beta }{{F^2}},\frac{r_a \sin \beta \sin {\rm \pi}\phi}{{F^2}}\right) . \end{equation}

Observe that the velocity gradients depend on two terms: a pressure gradient that is modified with the mean static pressure gradient, cf. (4.8), and then on a buoyancy term that scales with $\Delta \rho$. Note also that $\rho _1 = 1$, due to scaling with the density of fluid 1.

These expressions for the velocity gradients can be integrated further, outwards from the wall at $y=H$, where the velocity is zero, to give $v(y)$ and $w(y)$. A further integration over $[0,H]$ produces expressions for $(H\bar {v},H\bar {w})$. The point of this laborious exercise is to see that the integral expressions that multiply the modified pressure gradient terms and those that multiply the buoyancy terms, are identical in the two directions,

(4.14)\begin{equation} (H\bar{v},H\bar{w}) ={-}\mathcal{I}_1 \left(\frac{1}{r_a}\frac{\partial p}{\partial \phi}- \frac{\rho \,f_{\xi}}{r_a} , \frac{\partial p}{\partial \xi} + \frac{\rho \,f_{\phi} }{r_a} \right) + \mathcal{I}_2 \frac{\Delta \rho }{Hr_a}(\,{-}f_{\xi}, \,f_{\phi} ). \end{equation}

This is valid for any generalised Newtonian fluid. The expressions for $\mathcal {I}_1$ and $\mathcal {I}_2$ are defined as

(4.15)$$\begin{gather} \mathcal{I}_1= \left[ y_i \int^H_{y_i} \frac{\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y} + \int_0^{y_i} \int^{y_i}_y \frac{\tilde{y}}{\eta_2(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d}y} + \int_{y_i}^H \int^H_y \frac{\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} \right], \end{gather}$$

(4.16)$$\begin{gather}\mathcal{I}_2 = \left[ (H-y_i) \int_0^{y_i} \int^{y_i}_y \frac{\tilde{y}}{\eta_2(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} + y_i^2 \int^H_{y_i} \frac{H-\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y} + y_i \int_{y_i}^H \int^H_y \frac{H-\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} \right]. \end{gather}$$

These expressions depend on the constitutive laws for both fluids, as well as $y_i$ and $H$. For most generalised Newtonian fluids we will not be able to evaluate $\mathcal {I}_1$ and $\mathcal {I}_2$, except numerically.

Considering (4.14), in the absence of the buoyancy term, e.g. if $\Delta \rho = 0$, we can see that the mean flow follows the direction of the modified pressure gradient, as in Bittleston et al. (2002), and see (4.8). The mean density is again volumetrically averaged. With a density difference, however, the flow direction is influenced independently by the buoyancy vector: the second term on the right-hand side of (4.14). The form of (4.14) has the implication that the shear flow is a 1-D flow, in the direction of the gap-averaged velocity. Although we do not know this direction a priori, we may proceed assuming that this direction is known, solve the 1-D flow, then effectively calculate $\mathcal {I}_1$ and $\mathcal {I}_2$ to determine the flow direction from (4.14).

4.2.1. The 1-D shear flow

We assume that $(\bar {v},\bar {w})$ flows locally in direction $\boldsymbol {s} = (s_\phi,s_\xi )$. Taking the dot product of (4.5) and (4.6) with $\boldsymbol {s}$ we have the following 1-D momentum balance:

(4.17)\begin{equation} \frac{\partial }{\partial y} \tau_{sy} = \left\{ \begin{array}{@{}ll} \displaystyle{ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) - \dfrac{(H-y_i) \Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in [0,y_i],\\ \displaystyle{ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) + \dfrac{y_i \Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in (y_i,H]. \end{array} \right. \end{equation}

This type of buoyant two-layer flow model can generally be resolved numerically in an efficient way for two generalised Newtonian fluids. For example, see Zare, Roustaei & Frigaard (2017) for Bingham fluids, or Taghavi et al. (2009) for a different two-layer model for Herschel–Bulkley fluids. The problems are tractable due to various monotonicity results. Here, however, we focus on the case of Newtonian fluids for simplicity.

We first integrate to find the shear stress, which varies linearly in each layer,

(4.18)\begin{equation} \tau_{sy} = \left\{ \begin{array}{@{}ll} \displaystyle{\left[ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right]y - (H-y_i)y \dfrac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-f_{\xi}s_{\phi}) }, & y \in [0,y_i],\\ \displaystyle{ \left[ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] y - (H-y)y_i \dfrac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in (y_i,H] . \end{array} \right. \end{equation}

Note that we are not interested (yet) in the actual velocity solution $u_s(y)$, but in the areal flux $H\bar {u}_s$. This can be evaluated for Newtonian fluids by integrating by parts,

(4.19)\begin{align} H\bar{u}_s &={-} \int_0^{y_i} y \frac{\partial u_s}{\partial y}\,{{\rm d} y} - \int_{y_i}^H y \frac{\partial u_s}{\partial y}\,{{\rm d} y} ={-} \int_0^{y_i} y \frac{\tau_{sy}}{\eta_2}\,{{\rm d} y} - \int_{y_i}^H y \frac{\tau_{sy}}{\eta_1}\,{{\rm d} y} \nonumber\\ & ={-} \left( \frac{y_i^3}{3 \eta_2 } + \frac{H^3 - y_i^3}{3 \eta_1 } \right) \left[ \frac{\partial p}{\partial s} + \frac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] \nonumber\\ &\quad + \frac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi})\left[ \frac{(H-y_i)y_i^3}{3 \eta_2 } + \frac{y_i(H-y_i)^2(H+2 y_i)}{6 \eta_1} \right]. \end{align}

We may also evaluate $H\bar {u}_s$ from (4.14),

(4.20)\begin{align} H\bar{u}_s &= \boldsymbol{s} \boldsymbol{\cdot} (H\bar{v},H\bar{w}) \nonumber\\ & ={-} \mathcal{I}_1 \left[ \frac{\partial p}{\partial s} + \frac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] + \mathcal{I}_2 \frac{\Delta \rho }{Hr_a}(\,f_{\phi} s_{\xi} -\,f_{\xi} s_{\phi}) . \end{align}

Comparing these two expressions, we see that for a Newtonian fluid pair,

(4.21a,b)\begin{equation} \mathcal{I}_1 = \left(\frac{y_i^3}{3 \eta_2 } + \frac{H^3 - y_i^3}{3 \eta_1 } \right),\quad \mathcal{I}_2 = \left[ \frac{(H-y_i)y_i^3}{3 \eta_2 } + \frac{y_i(H-y_i)^2(H+2 y_i)}{6 \eta_1} \right] . \end{equation}

Power law and yield-stress fluid pairs will be dealt with in a subsequent study, requiring computation. Having solved for the 1-D shear flow, we may now go back to (4.14) and derive the modified 2-D model.

4.2.2. The D2DGA model

The gap averaged stream function is defined as ${\nabla }_a \varPsi = (2H\bar {w},-2H \bar {v})$, which can be substituted into (4.14). To derive the D2DGA model, we rearrange (4.14) to give two equations for the two pressure gradient components. We then cross-differentiate to eliminate the pressure,

(4.22)\begin{align} 0 &= \boldsymbol{\nabla}_a \boldsymbol{\cdot} \left[ \left( \frac{r_a}{2\mathcal{I}_1} \right){\nabla}_a \varPsi \right] + \boldsymbol{\nabla}_a \boldsymbol{\cdot}\left[ \left( \rho - \frac{\Delta \rho}{H}\frac{\mathcal{I}_2}{\mathcal{I}_1} \right) \boldsymbol{f} \right] \nonumber\\ & = \boldsymbol{\nabla}_a \boldsymbol{\cdot} [\boldsymbol{S} + \boldsymbol{b}] . \end{align}

We observe that (4.22) is an elliptic Poisson problem, of the same form as in Bittleston et al. (2002) and Maleki & Frigaard (2017). However, the assumption of a layered flow gives rise to significant differences.

Considering first $\boldsymbol {S}$, we see that for a Newtonian fluid,

(4.23)\begin{equation} \boldsymbol{S} = \left( \frac{r_a}{2\mathcal{I}_1} \right){\nabla}_a \varPsi = \frac{3r_a }{2H^3} \frac{{\nabla}_a \varPsi}{\displaystyle{\frac{\bar{c}^3}{\eta_2} + \frac{(1-\bar{c}^3)}{\eta_1}}},\quad \Rightarrow\quad \tau_w = \frac{3}{2H^2} \frac{|{\nabla}_a \varPsi |}{\displaystyle{\frac{\bar{c}^3}{\eta_2} + \frac{(1-\bar{c}^3)}{\eta_1}}} , \end{equation}

as opposed to $\boldsymbol {S} = [3 r_a\eta /2H^3] {\nabla }_a \varPsi$ for the 2DGA models. If the 2DGA approach were to yield the same $\boldsymbol {S}$, we see that a specific form of viscosity interpolation would be needed. In practice a linear interpolation has been used for 2DGA computations, with the assumption that the fluids are mixed.

Turning now to the buoyancy vector $\boldsymbol {b}$, we see that there are two parts in (4.22). Firstly, if we write $\rho = (1-\bar {c}) \rho _1 + \bar {c} \rho _2$, this helps us to see that the first part of $\boldsymbol {\nabla }_a \boldsymbol {\cdot} \boldsymbol {b}$ represents changes in the mean density in the direction $\boldsymbol {f}$. These buoyancy gradients are exactly those considered in the models of Bittleston et al. (2002) and Maleki & Frigaard (2017). The second part results specifically from the layered flow. If we write $\rho = \rho _2 + (1-\bar {c}) \Delta \rho$, we see that the buoyancy gradients from the mean density term have similar order of magnitude as those that result from the second part. In other words, the effect of the layered flow on buoyancy is significant. In more detail, we have

(4.24)\begin{equation} \boldsymbol{\nabla}_a \boldsymbol{\cdot} \boldsymbol{b} = \boldsymbol{\nabla}_a \boldsymbol{\cdot} \left[ \rho_2 + \Delta \rho \left( (1-\bar{c}) - \left[ \frac{m\bar{c}^3(1-\bar{c}) + \bar{c}(1-\bar{c})^2(0.5+\bar{c})}{m\bar{c}^3 + 1 - \bar{c}^3} \right] \right) \right] \boldsymbol{f} . \end{equation}

The gradients of $\boldsymbol {f}$ will typically be small as changes in the well orientation occur slowly.

With regard to transport of the fluid concentrations, we need to evaluate the fluxes. For Newtonian fluids this is straightforward, although algebraically long. The pressure gradients are eliminated using (4.14) and there is some manipulation using the expressions $\mathcal {I}_1$ and $\mathcal {I}_2$. This simplifies to

(4.25)$$\begin{gather} (q_\phi , q_\xi )=\frac{ r_a}{2} \left(- \frac{\partial \varPsi}{\partial \xi } , \frac{1}{r_a}\frac{\partial \varPsi}{\partial \phi } \right) \bar{c} \left[ \frac{m\bar{c}^2 + 1.5(1-\bar{c}^2)}{m\bar{c}^3 + 1-\bar{c}^3}\right] + \frac{\Delta \rho H^3}{6 \eta_2} \mathcal{I}_3(\bar{c},m)[{-}f_\xi,f_\phi ], \end{gather}$$

(4.26)$$\begin{gather}\mathcal{I}_3(\bar{c},m) = \frac{\bar{c}^2(1-\bar{c})^3[4m\bar{c} + 3(1-\bar{c})]}{2m[m\bar{c}^3 + 1 - \bar{c}^3 ]}. \end{gather}$$

We observe that the first term is very similar to the model of Maleki & Frigaard (2017), except that the transport term in this model is

(4.27)\begin{equation} \frac{r_a}{2} \left(- \frac{\partial \varPsi}{\partial \xi } , \frac{1}{r_a}\frac{\partial \varPsi}{\partial \phi } \right) \bar{c} . \end{equation}

We see that this is amplified by the concentration dependent flux term,

(4.28)\begin{equation} \frac{m\bar{c}^2 + 1.5(1-\bar{c}^2)}{m\bar{c}^3 + 1-\bar{c}^3}. \end{equation}

On inspection, we see that this amplification factor is the first term in (4.2), as derived by Yang & Yortsos (1997). Unlike these earlier studies, in the annulus the gap width and local mean velocity may now vary.

The second term in (4.25) is also related to (4.2). If for example we simplify to have a single component of $\boldsymbol {f}$, e.g. in a vertical well with $\beta = 0$, then we see that $f_\xi =0$. The $\phi$-component of $\boldsymbol {q}$ consists of only the first dispersive term of (4.25), whereas the $\xi$-component is augmented by buoyancy. The $\bar {c}$-dependence of $\mathcal {I}_3(\bar {c},m)$ mimics that of the second term in (4.2) and the expressions are identical if

(4.29)\begin{equation} U = \frac{3m \eta_2 {F^2}}{\Delta \rho H^3 r_a}. \end{equation}

If the annulus is uniform and concentric ($H=r_a=1$), we find that these terms are identical.

5. Results using the D2DGA model

Including the fluxes of (4.25) and solving (4.22) for the stream function constitutes the D2DGA model. With the corrected flux terms in the concentration equation, this takes the dispersion behaviour along the channel into account, but also modifies the effects of buoyancy and rheology on computing the velocity field. We start our comparisons by revisiting the computations of the total volume of displacing fluid in the annulus. We only show the mass conservation properties of cases 6 and 8, since earlier we found a distinct discrepancy between the 2DGA and 3-D results in terms of the breakthrough time, i.e. for $m>1$ the 3-D results have earlier breakthrough, see figure 4(b). Figure 11 shows the D2DGA model comparison with the 3-D results. The D2DGA model now also predicts earlier breakthrough times than the piston-like displacement and smaller volumetric efficiencies, very similar to the 3-D model. The inclusion of gap-scaled dispersion evidently leads to this improvement. Mass conservation was verified for other cases, with similar agreement; not shown here for brevity.

Figure 11. Verification of volume conservation for the D2DGA and 3-D models of cases 6 and 8. The black line represents the total annulus volume. The solid lines with other colours represent the theoretical value based on the constant flow rate. The D2DGA and 3-D results are marked by different symbols, as indicated.

We now turn to comparing the concentration fields at successive times. Based on the previous results, we focus first on the cases with relatively large discrepancy between 2DGA and 3-D models. Figure 12 shows the results of cases 7 and 3. Figure 12(a,c) display the 3-D results with the 2DGA results in the inset figures. Figure 12(a,c) show the results of the D2DGA model. We see a significant improvement of the D2DGA model compared with the 2DGA. Firstly, we observe clear concentration gradients, which represents the gap-scaled dispersion, which agree much better with the 3-D results. Secondly, we see instabilities near the displacement front on the narrow side for case 3 in both 3-D and D2DGA results. The underlying D2DGA model is still a Hele-Shaw-style model, except now the concentration dependence of $\boldsymbol {b}$ and $\boldsymbol {S}$ is significantly changed, as well as the transport terms. It seems that these changes allow for fingering-type instabilities to develop that are analogous to the gap-averaged 3-D instabilities.

Figure 12. Gap-averaged concentration colour maps of (a) case 7 computed by 3-D model; (b) case 7 computed by D2DGA model; (c) case 3 computed by 3-D model; (d) case 3 computed by D2DGA model. The 2DGA results are in the top left-hand corners as comparisons. Parameters: $e=0.2$, $Re=20$, $m=0.5$, $b=10$ (case 7); and $e=0.6$, $Re=20$, $m=0.2$, $b=10$ (case 3).

Although we cannot dismiss inertial effects, their role in these instabilities appears to be minor. We state this simply because there is no inertia in the D2DGA model, although of course the 3-D flows have inertial effects ($Re = 20$). Instabilities on the narrow side of the annulus have been observed experimentally (Tehrani, Bittleston & Long 1993), where the authors attribute them to density driven azimuthal currents. Here they appear to be fingering-like instabilities, which are usually due to buoyancy and mobility gradients between the fluids. Certainly azimuthal currents are present in both 2-D models, although they appear to be more amplified for the D2DGA model.

To pinpoint the origin of these instabilities is not easy. There is a direct influence on the D2DGA model from the fluxes in the concentration equation, i.e. direct dispersion. We can see the effect of this in the earlier time subpanels of figure 12(b,d), where the concentration spreads vertically. Since the annuli are eccentric, this vertical dispersion results in concentration gradients azimuthally, i.e. because the narrow side moves axially (and disperses) at a slower speed. We observe that it is not until significant azimuthal stratification emerges in figure 12(b,d), that instabilities emerge. This suggests that the onset comes from developing azimuthal flows, as previously postulated (Tehrani et al. 1993). Note that the elliptic equation for the stream function is driven by gradients in $\boldsymbol {b}$. In the case that the annulus is vertical these gradients arise from $\phi$-gradients of the concentration-dependent expression on the right-hand side of (4.24).

As discussed previously, the last two terms on the right-hand side of (4.24) have comparable size. We note that the first two terms are also present identically in the 2DGA model, but no such instabilities are observed for these cases. Arguably, however, since the 2DGA model does not disperse as much as the D2DGA model, in the concentration equation, the azimuthal gradients needed are simply slower to develop. Indeed Renteria & Frigaard (2020) do find azimuthal instabilities in the 2DGA model for cases when long narrow-side channels develop in horizontal wells. Thus, we believe that the instabilities result from a combination of both enhanced axial dispersion (creating the azimuthal concentration gradients) and the azimuthal gradients in $\boldsymbol {b}$, that generate secondary flows via $\varPsi$.

The only other change with the D2DGA model comes in the stream function equation, via $\boldsymbol {S}$. As we have discussed following (4.23), the change in $\boldsymbol {S}$ can be considered as equivalent to selection of a specific rheological closure for the mixture viscosity. In the 2DGA model we have used a linear interpolation for the rheology, but in its early development we had also explored different closure models and found minimal effect. Therefore, we suspect that the effects of $\boldsymbol {S}$ are not dominant.

The main discrepancy we find from cases 3 and 7 above is that the D2DGA model in fact gives a slightly more dispersive interface than that in 3-D, see also figure 11. The difference is more distinct on the wide side, where the 3-D results typically show a long thin dispersive spike profile, see figure 13 for a snapshot from case 3. At the front of the spike we observe that the tip becomes diffuse, effectively due to numerical resolution across the gap. The concentration is gap-averaged to give $\bar {c}(\phi,\hat {\xi },\hat {t})$ for the 3-D model. In figure 14 we have taken a slice along the annulus at $\phi = 0$ (wide side) and we plot $\bar {c}(0,\hat {\xi },\hat {t})$ for the three models. The curves are spaced at equal time intervals throughout the displacement, from which we can infer the front velocity. The 3-D and D2DGA models are qualitatively similar and quite different from the 2DGA model (see the inset in figure 14b), again confirming the relevance of gap-scale dispersion. At large $\bar {c}$ the front velocity is comparably slow, indicating that residual fluid at the wall is slow to remove. The front velocities then increase until attaining a plateau value of velocity (at around $\bar {c}=0.75$ in 3-D, and $\bar {c}=0.65$ in D2DGA). This front velocity is constant until low $\bar {c}$, when it again increases. This pattern results in the dispersive low concentration regions advancing ahead of a more substantial front, as we have seen in figure 12.

Figure 13. Spike penetration regime in wide channel gap ($\phi =0$) of case 3 modelled by 3-D.

Figure 14. Advective gap-averaged concentration in wide channel of case 3 (parameters: $e=0.6$; $Re=20$; $m=0.2$; $b=10$) modelled by (a) 3-D; (b) D2DGA (the 2DGA result is in the top right-hand corner).

For the analogous 1-D flows studied by Yang & Yortsos (1997) and Lajeunesse et al. (1999), i.e. (4.1) and (4.2), show a range of different qualitative behaviours. The profiles in figure 14 correspond to contact shocks (using the classification of Lajeunesse et al. (1999)), and indeed if we evaluate $m$ and $U$ for the parameters of case 3, this regime is indicated. The contact shock profiles have a region of small $\bar {c}$ that moves fast and ahead of the main front, i.e. disperses. Then there is a shock, over a range of intermediate concentrations, see figure 14(a). We adopt the terminology of Lajeunesse et al. in this paper mainly for simplicity.

Looking at the low $\bar {c}$ range of figure 14, the speed of the lowest $\bar {c}$ fronts appear faster in the 3-D case than the D2DGA, although these represent relatively little volume. We cannot expect an exact correspondence between the models. The advancing wide side front is also fed by azimuthal secondary flows in both models, the 3-D model still retains inertia and other stress gradients. Numerically, for the D2DGA we resolve a model that is derived by integrating a ‘sharp interface’ on the gap-scale, to give closure expressions for $\boldsymbol {b}$ and $\boldsymbol {S}$. This is compared with the 3-D model where we gap-average the result, e.g. figure 13, which becomes diffuse as it advances and thins to the thickness of mesh used across the gap.

As well as contact shocks (spike), Lajeunesse et al. (1999) also identify characteristic behaviours termed no shock (fully dispersive) and frontal shock (shock), found in well-defined regions of the ($m$–$U$) plane. We can find also these other characteristic behaviours in our results. Figure 15(a) shows the concentration $\bar {c}(1,\hat {\xi },\hat {t})$, along the narrow side of case 7 ($m=0.5$), which has no shock (dispersive). Figure 15(b) shows $\bar {c}(1,\hat {\xi },\hat {t})$ for case 6 ($m=5$), which has a frontal shock. These regime shifts are quite sensitive to the viscosity ratio $m$ at fixed $U$ (i.e. $b$). For other parameters we can also find contact shocks on the narrow side, as well as experimentally (Etrati & Frigaard 2019). We have focused on wide and narrow side here, where the 2-D models have velocity only along the annulus (due to the imposition of symmetry).

Figure 15. Advective gap-averaged concentration profiles along the narrow side of the annulus, modelled by D2DGA for (a) case 7 (parameters: $e=0.2$; $Re=20$; $m=0.5$; $b=10$); (b) case 6 (parameters: $e=0.4$; $Re=20$; $m=5.0$; $b=100$).

Finally we present the concentration colour map of cases where there is little instability. Case 6 in figure 16 displays a typical frontal shock behaviour. The D2DGA result agrees remarkably with the averaged 3-D concentrations. Both models predict a much less dispersive interface compared with the previous cases shown. The flow exhibits a consistent degree of dispersion as the front propagates steadily along the annulus, with a modest gradient from wide side to narrow side, which is not captured at all in the 2DGA model.

Figure 16. Gap-averaged concentration colour maps of case 6 (parameters: $e=0.4$; $Re=20$; $m=5.0$; $b=100$) computed by (a) 3-D model; (b) D2DGA model. The 2DGA result is in the top left-hand corner as a comparison.

Case 4 is a buoyancy-dominated flow that showed similar steady front speeds between 2DGA and 3-D results (see figure 5). It can be seen from figure 17(b) that the D2DGA result also exhibits a stable/steady front with little dispersion on the wide side. A similar degree of narrow side dispersion in predicted to the 3-D result. In both cases 6 and 4, we can discern some faint wave-like patterns in the averaged 3-D solution at later times. As with our earlier detailed examination of case 7, these waves arise from a loss of azimuthal symmetry behind the front. This asymmetry is prevented in the D2DGA model as implemented, since we have imposed a symmetry condition on $\varPsi$ at the wide and narrow sides. However, relaxing this constraint (replacing with periodicity) would allow for the asymmetry and is a small change computationally. The 2.5-D model of Tardy (2018) does appear well able to represent azimuthal asymmetry and at its foundation this model is also non-inertial. We cannot, however, discount that there may also be an inertial component to the symmetry breaking in the 3-D case.

Figure 17. Gap-averaged concentration colour maps of case 4 (parameters: $e=0.6$; $Re=20$; $m=0.2$; $b=1000$) computed by (a) 3-D model; (b) D2DGA model. The 2DGA result is in the top left-hand corner as comparison.

Taking into account that the dispersion behaviour appears to improve the 2DGA significantly. In particular, we may now use this model to predict the breakthrough time as well as the displacement efficiency more accurately than before, and also less time-consuming compared with 3-D models. This allows us to explore the effects of different dimensionless parameters in a more intuitive way through these two factors. We define a dimensionless breakthrough time as $t_{br}=\hat {t}_{br}\times\hat {w}_{0}/\hat {L}_{annu}$, and the displacement efficiency as $\eta _E$ which is calculated by dividing the fraction of the annulus displaced at a specific time ($1.2\times\hat {L}_{annu}/\hat {w}_{0}$) by the annulus volume.

Table 3 shows $t_{br}$ and $\eta _E$ as calculated by the D2DGA and 3-D models. The buoyancy number is the dominant parameter, as demonstrated by a later breakthrough time and higher efficiency with increasing $b$. We verify again via the poor efficiency of case 1 that $b <0$ should be strictly avoided. In addition, it seems that we find an earlier breakthrough time predicted by the 3-D model, when $b$ is small ($b \leq 10$), which corresponds to the long thin spike profile observed in figure 12(a,c). Note that the displacement efficiency predicted by both models are pretty similar for these cases. The effect of eccentricity is not very apparent through comparison between cases 2, 5 and 9, since the buoyancy number is large enough to compete against the geometry influence. However, it is worth noting that a small residual layer on the narrow side can result in extremely ineffective cementing, from the perspective of well integrity/leakage, even when $\eta _E$ is large. Thus, we should always combine $\eta _E$ with examination of the fluid placement to better determine the quality of a displacement. We can also use the D2DGA model to calculate for example a narrow side efficiency $\eta _N$, which focuses only on the most problematic section of the annulus and often gives a better representation of the efficacy of the displacement. There are various ways to do this (Guillot, Desroches & Frigaard 2007; Maleki & Frigaard 2018; Jung & Frigaard 2022).

Table 3. Breakthrough time and displacement efficiency of cases 1–10.

The effect of the viscosity ratio becomes apparent via comparison between cases 5 and 6, and cases 7 and 8, computed by the D2DGA model. Increasing the viscosity ratio causes more dispersion along the annulus and leads to lower efficiency, especially for cases with smaller buoyancy number ($b=10$). Hence, using fluid pairs with smaller viscosity ratio (displaced/displacing) is beneficial for complete and efficient displacement, as is recommended in industrial practice (Nelson & Guillot 2006). Interestingly, focusing on cases 5 and 6, we found that the viscosity ratio does not much affect the efficiency predicted by the 3-D model, unlike that of the D2DGA model. We suppose that 3-D effects, e.g. flow in radial direction, could compensate the severe dispersion caused by large viscosity ratio, specifically when the buoyancy force is strong. Finally, it seems that the Reynolds number also plays a less important role regarding the scale of dispersion and the displacement efficiency, (compare cases 5, 9 and 10). This observation is confined to the laminar regime flows studied. If turbulent displacements are considered then the dispersive picture changes substantially (Maleki & Frigaard 2017, 2018, 2019).

6. Summary

This paper has studied primary cementing displacement flows in vertical eccentric annuli. These flows are challenging to compute, due to the excessively long length scales of typical wells (${\sim} 10^3$ m), compared with the annular gap (${\sim} 10^{-2}$ m). Ensuring a mesh size needed for accuracy, both in representing the interface and on the scale of the annular gap, while also integrating in time long enough for the displacement front to travel the length of the well, makes 3-D non-Newtonian computations unfeasible at present. Nevertheless, it is important to understand the effectiveness of the full displacement process for engineering design purposes. By design we mean selection of the annulus eccentricity, the fluid densities and rheologies, the volumes and flow rates pumped. Each of these may only be varied within operational constraints. Because of this need, 2DGA models have become popular over the past 20 years.

In this paper we have used 3-D transient computations in annuli of length equal to ${\sim}10^3$ annular gap widths, to represent the full complexity of the displacement flow, and compare these with 2DGA results. While 2DGA models can predict many qualitative phenomena well, (e.g. steady/unsteady, stable/unstable, residual narrow side channel or not), there are also quantitative discrepancies in most cases. Computed 3-D simulations predict fluid concentrations that are clearly dispersive and that result in lower displacement efficiencies than 2DGA flows.

To account for gap-scale dispersion in the 2DGA framework, we rederive the model, starting from the reduced shear flow equations. We assume that the displacing fluid can symmetrically displace down the centre of the local annular gap, instead of assuming a fully mixed concentration across the gap (Bittleston et al. 2002; Maleki & Frigaard 2017). This leads to a buoyant two-layer flow model for the 1-D shear flow. The D2DGA model is obtained after solving the shear flow equations and eliminating the pressure components. This model derivation is presented for any two generalised Newtonian fluids. However, for many of these a numerical solution on the gap-scale would be needed. On the other hand, for the two Newtonian fluids that are used in our examples, the model can be found analytically.

The D2DGA equation for the stream function is an elliptic Poisson problem, just as in the 2DGA model. However, there are differences in closure expressions for both the wall shear stress and buoyancy. With regard to the mass transport equation, we obtain a newly defined gap-averaged flux that can approximate the actual dispersive behaviour. An additional term representing buoyancy effects in the layered flow also arises, compared with the 2DGA model. As the underlying mathematical structure of the 2DGA and D2DGA models is the same, computational times are similar, with some additional cost from evaluating the new closure functions.

Of course, the D2DGA model is not the first study of dispersive effects in duct-like displacements. In the classical study of Saffman & Taylor (1958) the authors explore viscosity (mobility) ratio effects on instabilities in immiscible Hele-Shaw displacements: viscous fingering. Much of this paper concerns viscous fingering in the case in which there is a single fluid present in the through-thickness direction. However, the Appendix contains a modified analysis of the situation with incomplete displacement, i.e. the displacing fluid has a multilayer structure. As here, this renders the mean velocity concentration dependent. Secondly, multiphase porous media flow models applied in oil recovery often employ flux functions that describe the wetting of different phases, dating from Buckley & Leverett (1942). These expressions are either empirically derived or rely on simplified models, such as the Hele-Shaw cell or tube bundles, in order to develop mathematical closure expressions. The wetting closures are then applied to rather complex actual pore geometries. Thus, our approach is certainly related to these larger application areas.

Perhaps we are fortunate in that the narrow eccentric annulus of primary cementing forms a prototype that has systematic and understandable geometry variation (effectively the permeability), allowing for the modelling approach to be effective. As displacement flows in cementing are also (predominantly) miscible and large Péclet number, our extension of Yang & Yortsos (1997) and Lajeunesse et al. (1999) seems the appropriate modelling approach. For uniform, vertical concentric annuli, the new closure expressions we derive are identical with those of Yang & Yortsos (1997) and Lajeunesse et al. (1999), who considered high Péclet number miscible Newtonian displacement flows in long vertical channels.

Through detailed comparisons between the D2DGA model and the 3-D model, we find significant improvement of the D2DGA model in capturing gap-scale dispersion. Moreover, some instabilities are predicted by D2DGA model, similar to those in the 3-D model. The D2DGA model remains computationally efficient and appears able to represent the different scales of dispersion, consistent with (averaging of) 3-D simulation results. The one deficiency is in modelling density-unstable displacement flows. Although the D2DGA prediction may be an improvement, quantitative prediction remains poor, as might be expected.

Finally, the effects of the dimensionless parameters have been explored, in terms of breakthrough times and displacement efficiency. The buoyancy number is the key factor in vertical primary cementing, once the eccentricity is fixed. Negative buoyancy number should be strictly avoided to prevent highly ineffective cementing and compromised well integrity. A large enough buoyancy number can be effective for Newtonian fluids, even if the annulus is highly eccentric. It helps with displacing the fluids on the narrow side and reducing residual layers caused by eccentricity. Small viscosity ratio (displaced/displacing) is also advisable. In the D2DGA model this directly decreases the dispersion between displacing and displaced fluids, essentially by acting on the drainage layers adjacent to the walls. Especially when the buoyancy number is small, the viscosity ratio affects dispersion more significantly, i.e. large viscosity ratio $m$ will lead to higher dispersion and less efficiency. Eccentricity affects the secondary flow and thus contributes to uniform/non-uniform distribution in the azimuthal direction. Large eccentricity is not beneficial for complete and efficient displacement. Here we only studied Newtonian fluids, but for yield stress fluids it is also well known that these fluids can become stuck on the narrow side for large $e$, (McLean et al. 1967; Couturier et al. 1990). As for the effect of Reynolds number, we did not find that it plays an essential role and especially with large $b$. Although operationally speaking, this simply affirms an old story (Couturier et al. 1990; Frigaard & Pelipenko 2003; Pelipenko & Frigaard 2004c; Nelson & Guillot 2006), this paper gives far greater detail. The real novelty is in being able to access predictions of the dispersive behaviour that agree well with the more expensive 3-D models.

The next step in our study is to consider how effective the D2DGA model is in predicting experimental results from our flow loop. Thereafter, the model will be applied to inclined and horizontal annuli, which can be achieved easily by changing the direction of the buoyancy force ($\beta \neq 0$). Here we already have the recent studies of Renteria & Frigaard (2020) and Sarmadi et al. (2021) to compare with. The next step is to develop the model targeting non-Newtonian fluids, as are more common in field situations with laminar flows. Although a more complicated computation is needed to resolve the gap-scale flow, e.g. for Herschel–Bulkley fluids, the basic principle is applicable as outlined in this paper. Another interesting question is whether the new buoyancy, wall shear stress and flux closures can be analysed in a similar way to Pelipenko & Frigaard (2004c) and Maleki & Frigaard (2020), to rapidly classify cementing displacements in a way that is suitable for machine learning approaches.