A Numerical Study of Density-Unstable Reverse Circulation Displacement for Primary Cementing

The removal or displacement of one fluid by another is common in the petroleum industry and in chemical and process industries, as well as in nature. The cementing operation is an important example from petroleum industry, where the annular space between two casing strings of different diameters, or between a casing string and the rock formation, is to be displaced to a cement slurry. The annulus is initially filled with a drilling fluid, and a sequence of wash and spacer fluids is typically injected before the cement slurry. An important function of the wash and spacer fluids is to completely remove the drilling fluid in order to avoid inter-mixing with the cement slurry. Once pumped in place, the slurry will harden into a solid cement sheath, providing mechanical support for the cemented casing string and preventing uncontrolled migration of formation fluids behind the casing [1].

The conventional method of primary cementing involves injecting spacer fluids and cement slurry into the casing string that is to be cemented. An illustration of conventional circulation displacement is provided in Fig. 1(a), with the left schematic indicating fluid locations during injection and the right schematic showing intended fluid locations at the end of the operation. As shown in the left schematic, the injected fluids first displace the drilling fluid from the inside of the casing. The injected fluids exit the casing string at the bottom and continues up towards the surface within the annulus behind casing, where displacement of drilling fluid continues. The injected fluids will normally have different mass densities, where spacer fluids tend to have intermediate density between the drilling and the cement slurry, which is the denser fluid. Washing fluids are often low-density brines or oils that are pumped in the turbulent regime. The hierarchy of fluids when pumped down the well inside the casing is therefore density-unstable, with denser fluids above the lighter wash or drilling fluid. Inter-mixing between fluids can be avoided inside the casing by utilizing at least one elastic plug, or wiper plug, that provides physical separation of the fluids and mechanical washing of the inner wall of the casing [1]. Once in the annulus, the hierarchy is again mechanically stable, at least with respect to the cement slurry which is now below the other fluids. The pumping stops when the top plug lands on top of the bottom plug, indicated in the right schematic in Fig. 1(a). Now, the slurry is allowed to harden into a solid sheath that should provide zonal isolation behind the casing [2]. Some wells develop sustained casing pressure (SCP) or surface casing vent flows over time [3], which could be due to, for example, mud contamination [4], residual wall layers [2,5], annular gas migration [6], or thermal and mechanical loads during construction or production [7–10].

A combination of laboratory experiments, numerical simulations, and yard tests has advanced the study and knowledge of conventional circulation cementing and mud displacements to a relatively mature stage (see, e.g., the review by Frigaard et al. [11] and references therein). Design rules for effective laminar flow displacements have also been proposed, recommending minimum density and friction pressure hierarchies among the fluids to enable stable and steady displacements [12–14]. Although conventional circulation cementing is the more traditional and common approach to primary cementing, it can in some cases be beneficial or even necessary to reverse circulate cement in place. Toward the end of the reverse cementing operation, the narrow, eccentric annulus is filled with a viscous cement slurry that exerts a relatively high friction pressure [15–17]. The casing is filled normally with less viscous drilling fluid and spacer fluids, resulting in a smaller friction pressure. Consequently, the circulating pressure toward the bottom of the hole during reverse circulation cementing can be reduced compared to conventional displacements. This can enable placement of more viscous and denser cement slurries compared to what is possible by conventional circulation. Reducing the need for retarders in the cement slurry is another perceived benefit of reverse circulation cementing [18–21]. In particular, reverse circulation has been a preferred strategy for cementing certain high-temperature geothermal wells in Iceland, where long cemented casing sections are required in order to withstand extreme thermal and mechanical loads [22]. Cementing by reverse circulation involves injecting the cement slurry directly into the annulus from the surface and displacing drilling fluid and other fluids downward, as illustrated in Figs. 1(b) and 1(c). Access to the annulus and cementing by reverse circulation are normally possible in onshore wells and platform wells, and solutions for offshore wells with subsea wellheads have been proposed [23], making this a relevant placement option for an increasing number of wells. Occasionally, reverse circulation cementing may also become a necessary second stage of a cementing operation, for instance, if a weak zone prevents completing the cementing job as a one stage operation. This is illustrated in Fig. 1(c). Here, a first stage conventional circulation operation should ensure cementing of the casing up to the loss zone. The goal of the reverse circulation second stage is to cement the casing interval above the loss zone by bullheading cement slurry down the annulus from the surface.

Where conventional circulation cementing places the cement slurry below drilling fluid and spacer fluids in the annulus, and thereby normally leads to density-stable configurations, reverse circulation places the cement slurry above the normally lighter fluids in the annulus. This mechanically unstable configuration is likely more prone to develop instabilities and inter-mixing of fluids during cementing compared to conventional circulation operations. The goal of this study is to start exploring density-unstable annular displacements, using insights from a North Sea production well where the top part of the production casing was originally cemented in 1985 following a procedure similar to that shown in Fig. 1(c). The annular cement behind the top part of the production casing in this well has been the topic of active study in recent years, see, e.g., Refs. [24–27]. The recent research activities stem from the 2018 well abandonment operation and the retrieval of two sandwich joints consisting of production casing, annular cement, and intermediate casing as part of this operation. The two retrieved joints have been re-logged at the surface using modern ultrasonic logging tools [24,25], and they have been the subject of pressure and seepage tests for annular permeability characterization [25,26]. Interestingly, the upper part of the production casing was cemented by reverse circulating (or bullheading) cement slurry down the annulus from the surface. The bullheading operation was performed as the second stage cementing job, following the first stage conventional circulation cementing of the bottom part of the production casing which covered the casing shoe and up to a weak formation layer. Although the well did develop a history of sustained casing pressure (SCP) and evidence of fluid migration behind the production casing, visual inspection of the annular cement in the two sandwich joints showed complete cement coverage around the entire annulus and relatively few and small inclusions of foreign material in the hardened cement. The exposed annular cross sections of the two sandwich joints are shown in Fig. 2, from Ref. [25]. One of the joints, referred to as the transition joint, contains the top-of-cement in the middle of the joint. The material toward the top of this section is settled barite, probably from the fluid pumped behind the cement slurry. Both surface re-logging and the ensuing permeability characterization concluded that the top-of-cement is relatively sharp and well-defined, indicating minor mixing and contamination at this location. Subsequent pressure and seepage testing of the joints has indicated permeabilities corresponding to effective microannuli smaller than 30 microns, which is safely within the interval of representative permeabilities for wells that exhibit SCP [28]. This suggests that the reverse circulation cement operation was relatively successful, even if it failed to provide complete zonal isolation.

The density-unstable fluid configuration associated with reverse circulation cementing is connected to the classic Rayleigh–Taylor instability [29,30], albeit in the highly confined annular geometry behind a casing string and in the presence of an imposed axial flow. Effects of confinement on the Rayleigh–Taylor instability were studied by Debacq et al. who performed buoyant exchange flow experiments in a vertical, slender pipe geometry using two miscible fluids [31]. The study showed that cross-sectional averaged concentration profiles at different times $t^$ since the beginning of the experiment approached a single, collapsed curve when the position $x^$ was scaled according to the similarity variable $x^/t^$⁠, equivalent to the similarity coordinate in microscopic diffusion processes. Here, $x^$ is an axial position relative to the initial, flat interface between the fluids, and the shape of the collapsed curve was determined by a parameter referred to as the “macroscopic” diffusion coefficient of the buoyant displacement. Further experimental work connected to density-unstable exchange flows in pipes at inclinations ranging from vertical to near-horizontal was performed by Séon et al. who further studied macroscopic diffusion [32], front velocity [33,34], and front dynamics [35] within exchange flows in pipes. Buoyant displacements were classified as either macroscopically diffusive, inertial, transitional, or viscous, depending on the dominating balances of inertia, viscous, and bouyant stresses in the experiments, as well as pipe inclination [35].

In a perfectly vertical pipe, effective transverse mixing will limit the front velocity of the mixing zone. As inclination from vertical increases, the dense fluid will increasingly tend toward the low side of the pipe, which in turn increases downward and upward front velocities of the mixing zone [35]. These exchange flows are inherently unstable, resulting in a growing mixing zone or stratification of the fluids. An imposed flow velocity may counteract the tendency for the lighter fluid to flow upward, and in this sense has a stabilizing effect on the density-unstable configuration. Taghavi et al. [36,37] and Alba et al. [38] studied density-unstable displacements with downward imposed flow in pipes and in channels, using both experimental and computational approaches. In Ref. [38], a comprehensive displacement regime classification based on a significant experimental database was presented, categorizing displacements as either backflow or instantaneous depending on whether traces of lighter fluid appeared above the initial interface position or not. Displacements were further categorized as macroscopically diffusive, viscous, or inertial [38]. More recently, Amiri et al. investigated miscible iso-viscous and density-unstable downward displacements in a slender, vertical pipe, and categorized displacements as stable or inertial and identified conditions for front detachments and decaying helical waves [39]. A 1D lubrication model was also developed to support displacement categorization, assuming a downward core-annular flow with the dense fluid occupying the central core region. Model predictions compared favorably to observations, especially when viscous stresses balanced buoyancy stresses [39]. Recent work has also studied effects of non-Newtonian viscosities, such as the study of density-unstable displacements in inclined pipes by Alba and Frigaard [40], where a denser Newtonian fluid displaced a yield stress and shear thinning fluid downward. As most well construction fluids are multi-component fluids that normally exhibit both yield stress and shear thinning properties, as well as develop a significant gel strength under static conditions [41], rheological effects will impact displacements in wells. Residual wall layers of unyielded, static fluid can be particularly difficult to remove with only the shearing action of displacing fluid [40].

Numerical experiments and simulations of density-unstable exchange and displacement flows complement physical experiments and can help explore details that are not immediately accessible in experiments. Buoyancy-driven mixing in two dimensions and three dimensions were compared numerically by Hallez and Magnaudet [42], using direct numerical simulations that were benchmarked with the experimental data of Séon et al. [34]. In line with the remark above, the simulations were considered to provide an improved connection between the evolution of the concentration profiles observed in experiments and the internal shear field in the fluids [42]. Comparing results obtained 3D geometries and the equivalent, but more confined 2D geometry suggested significant differences, such as in the front velocities, which were larger (and in better agreement with experiments) in 3D simulations [42]. Subsequent studies of exchange flows include a detailed numerical analysis of secondary velocity components in the mean flow [43] and an integrated experimental–numerical analysis of momentum transport and turbulent structure [44]. More recently, Hallez and Magnaudet presented a comprehensive analysis of density-unstable exchange flow in a near-vertical pipe, using the experimental configuration in Ref. [44] as the basis for the direct numerical simulation [45]. The study provided further insights into turbulence production and transport, and effects of the pipe walls on, for example, energy redistribution. Finally, Sebilleau et al. explored the performance of different numerical solvers for simulating such density-unstable exchange flows [46].

While the above studies have explored density-unstable exchange and displacement flows in pipes and channels, the displacement mechanics of annular reverse circulation cementing have not been studied in detail previously. The purpose of this paper is to utilize the available data from the original 1985 reverse circulation cementing operation discussed above to study the density-unstable displacement process in the well using fully 3D computational fluid dynamics (CFD) simulations. The outline of the paper is as follows. The governing equations and relevant dimensionless numbers are introduced in the next section, along with the parameters used in the numerical simulations. In the subsequent section, results are presented in terms of gap-averaged and azimuthally averaged concentration profiles for the two inclinations. The results are compared with a simplified, 1D displacement model for plane channel displacements. Finally, we discuss the impact of the results for interpreting the sandwich joint cement quality before concluding. Details of the 1D model derivation are provided in the  Appendix.

This study is focused on density-unstable annular displacements involving two iso-viscous Newtonian fluids that are considered incompressible. The two fluids are initially at rest and separated by a sharp interface that is parallel to the cross-sectional plane of the annulus with the denser fluid above the lighter fluid. A constant flux of displacing fluid is enforced through the annular cross section at the inlet, while the fluids flow freely out of the lower outlet cross section.

In the simulations that follow, the fluid densities and viscosities, and the dimensions of the annulus remain fixed while the imposed bulk velocity is varied. This motivates using the imposed velocity as a scale for the velocity. Other velocity scales can be identified for the problem at hand by balancing viscous and buoyant stresses, or inertial and buoyant stresses. The former would suggest a velocity scale $U^v*=Atg^D^h*2/ν^$⁠, with $ν^=μ^/ρ¯^$ when the flow is dominated by viscous and buoyant stresses. The latter, balancing inertial stresses $ρ¯^U^i*2$ by buoyant stresses $Δρ^g^D^h*$⁠, suggests $U^i*=Atg^D^h*$ as characteristic velocity scale in this case.

The basis and initial motivation for this study is the reverse circulation cementing operation of the production casing in the North Sea production well discussed above. The geometry of the annulus is given by the outer diameter of the 9 5/8-in. production casing (corresponding to a radius of 122.2 mm) and the inner diameter of a 13 3/8-in. intermediate casing (corresponding to a radius of 156.8 mm). This produces a diameter ratio of approximately 0.78 and a relatively narrow annulus. The simulations presented in this study are all from a model geometry (illustrated in Fig. 3) that is 8 m long in the axial direction, and with the initial interface situated 2 m below the inlet at the top. The total axial length corresponds to approximately 116 hydraulic diameters.

For the well in question, the upper part was cemented by reverse circulating cement in place as a second stage of the cementing operation. Consequently, near-vertical inclinations are considered most relevant. In the following, a strictly vertical annulus and an annulus inclined by 10 $deg$ from the vertical will be studied. Furthermore, at both inclinations, the study will focus exclusively on concentric configurations. While casings are in practice never perfectly centralized, even in perfectly vertical wells, we choose to limit the current study to focus on the effect of inclination. As will be suggested below, eccentricity and inclination are anticipated to have similar effects on the displacement flow: in the case of inclination, the component of gravity that is transverse to the axis of the annulus will accelerate dense, displacing fluid toward the low side of the annulus. In the case of a vertical eccentric annulus, the variable gap width will now direct the dense fluid toward the wider sector of the annulus. In both cases, the preferred direction of flow is anticipated to elongate and destabilize the interface between the fluids.

For the two fluids involved in this study, their mass densities are taken as 1820 kg/m³ and 1900 kg/m³ for the displaced and displacing fluids, respectively. Both densities have been determined on the basis of the original 1985 cementing job reports. Unfortunately, little is known with certainty about the specific viscosities or injection rates used during cementing. The composition of the slurry was very simple with no additives specified. Consequently, and inspired by measurements on another simple well cement slurry [47], a constant viscosity of 0.1 Pa s is taken for the displacing fluid in the simulations. In the absence of any data for the spacer fluid or drilling mud used in the original cementing operation, the viscosity of the displaced fluid is also set to 0.1 Pa s for simplicity.

For the injection rates to be used in the simulations, a constant rate of 800 l/min is taken as the base case, again inferred from the original job report. This corresponds to an imposed bulk velocity of approximately 0.44 m/s in the annulus. In addition to studying this relatively high imposed velocity, the effects of lower imposed velocities on the displacement will also be studied by using bulk velocities of 0.1 m/s, 0.0437 m/s, and 0.01 m/s. The velocities correspond to flowrates of 182 l/min, 80 l/min, and 18 l/min in this annular geometry. A summary of the model parameter values is provided in Table 1, where the dimensionless number χ = 2Re/Fr² has been introduced. This number reflects the ratio of buoyant stresses to viscous stresses, and occur in the simplified 1D displacement model discussed in Sec. 3.3. Finally, one may compare the dimensionless numbers listed in Table 1 with the regime classification for density-unstable pipe flows in Ref. [38]. The combination $Re/Fr=55$ and the four values of Fr place these displacements in either a viscous regime (Fr = 1.2, 5.1) or an inertial regime with backflow (Fr = 0.1, 0.5). A comparison with this regime classification will be performed below, when analyzing results of the numerical simulations.

The governing equations are solved using the open-source computational fluid dynamics platform openfoam version 1912 and its volume-of-fluid solver interFoam, without surface tension between the fluids. Second-order discretization schemes are used for the spatial terms, while the so-called Crank–Nicolson scheme is used for the time discretization. The latter is a blend between the second-order Crank–Nicolson scheme and the first-order Euler scheme, implemented in openfoam with an option to enhance robustness by tuning a numerical factor, where 1 corresponds to pure Crank–Nicolson and 0 is pure Euler. A value of 0.9 typically ensures sufficient stability and will be used here.

For the boundaries, as shown in Fig. 3, velocity Dirichlet conditions are prescribed at the walls and at the inlet (no slip and uniform bulk inlet velocity, respectively) and a Neumann outflow condition at the outlet. For the pressure, Neumann inlet and Dirichlet outlet conditions are used.

While the imposed velocities correspond to a range of Reynolds numbers that indicate laminar flow, the buoyancy will have a significant destabilizing effect. To limit the scope of the current study, no turbulence model will be applied; however, a mesh with 30, 120, and 1024 cells in the radial, azimuthal, axial directions, respectively, as shown in Fig. 3(b), is found to yield sufficiently high accuracy for the displacement front velocities when compared to finer meshes. Previous validation of the numerical setup and solver were performed by comparing simulation results to conventional circulation displacement experiments in annuli at near-horizontal [48] and near-vertical inclinations [14]. Furthermore, as shown below, predictions of front velocity from the full 3D simulations agree favorably with results from a simpler 1D lubrication scaling model within the relevant range of viscous, density-unstable displacements (see Sec. 3.3).

While the numerical simulations are fully 3D, the result presentation that follows is based on post-processing of the results by averaging the local fluid concentration field in one or in two spatial directions. The main motivation for spatial averaging is to visualize the dominating trends in the simulations, which can be less clear in a full 3D rendering of the simulation results. To be specific, the fluid concentrations are averaged either radially across the gap between the walls of the annulus or azimuthally around the annulus. This step generates 2D diagrams, showing the averaged fluid concentration as a function of axial and azimuthal position when averaging radially. In the case of averaging in the azimuthal direction, 2D diagrams of averaged concentration as a function of axial and radial positions result. These separate averaging procedures are illustrated in Fig. 4. The gap-average (1) in Fig. 4 is equivalent to “unwrapping” the annulus and radially averaging the concentration field across the gap between the inner and outer walls. The azimuthally averaged concentration (2) is also equivalent to an “unwrapping” of the annulus, but the average is now taken in the azimuthal direction. The purpose of this averaging method is to resolve the average location of the fluids within the radial gap.

An example of these diagrams is provided in Fig. 5, where flow is from top toward bottom. The two pairs of diagrams shown here are taken at two different points in time, for a concentric and inclined annulus. The diagrams are labeled either (1) or (2), corresponding to the two different averaging directions shown in Fig. 4. Furthermore, in diagrams corresponding to (1), the vertical direction along the diagram corresponds to axial position while the horizontal direction corresponds to azimuthal position. As indicated by the top-most annotations, the left corner of the diagrams corresponds to the lower side, or bottom (B), of the annulus as seen in an inclined wellbore. The middle corresponds to the top (T) and the right corner is again at the bottom (B). For the diagrams labeled as (2), the vertical direction again corresponds to axial position. Now the horizontal direction in the diagram corresponds to radial position, from the inner wall of the annulus at the left edge (⁠$R^i$⁠) to the outer wall at the right edge of the diagram (⁠$R^o$⁠).

Finally, a second type of diagrams referred to as spatio-temporal diagrams will also be presented below. In this case, the fluid concentration is averaged in both directions, i.e., both radially and azimuthally, to generate a full cross-sectional average at each axial position at each time-step of the simulation. Plotting the averaged concentrations as a function of time provides quantitative information about the extent and evolution of the fluid mixing zone, as well as qualitative features within the mixing zone. Furthermore, the full cross-sectional average will also be used to identify possible macroscopic diffusion behavior, i.e., to assess whether the averaged concentration profiles at different times are similar when plotted as a function of the similarity variable $(x^−x^0−U^t^)/t^$⁠, where $x^0$ denotes the initial position of the horizontal fluid interface, and $U^$ the magnitude of the imposed bulk velocity.

The diagrams in Fig. 6 show the averaged fluid concentration profiles at selected times since start of the simulation for the four imposed bulk velocities in a vertical, concentric annulus. Two diagrams are presented at each time with the left showing a radially averaged concentration map (corresponding to the gap-average (1) in Fig. 4) and the right being the corresponding azimuthal-average (corresponding to (2) in Fig. 4).

Starting with the right-most diagrams corresponding to simulations where Fr = 5.1 or $U^=0.44$ m/s, a steady downward displacement is observed. The two diagrams show little or no variation in fluid concentrations in the azimuthal direction. The right diagram at Fr = 5.1 suggests the displacement up to this point in time and is largely driven by the velocity profile, i.e., being close to simple advection of the fluids with little influence of buoyant stresses.

The destabilizing effect of buoyancy becomes more pronounced at lower values of Fr, corresponding to lower imposed bulk flowrates. Already at Fr = 1.2, instabilities are present along the mixing zone. This is visible in terms of variations in the azimuthal direction, and also in that the azimuthally averaged diagram shows a less distinct interface compared to Fr = 5.1. This suggests more fluctuations in the interface shape in the circumferential direction. There is a clear tendency for the denser fluid to propagate downward through the center of the radial gap at both Fr = 1.2 and Fr = 5.1.

Turning next to the more density-unstable results at the lowest imposed velocities, corresponding to Fr = 0.1 and Fr = 0.5, variations in radial and azimuthal directions increase further. In both simulations, there is backflow of the lighter fluid, i.e., some of the lighter moves counter direction to the imposed flow and ends up above the initial interface position. Backflow is evident at the lowest bulk velocity, Fr = 0.1, where a significant volume fraction of lighter fluid travels upward. At Fr = 0.5, this trend is clearly diminished, although the fluid configuration at $t^=30$ s in the figure reveals a slight migration of lighter fluid above the location of the initial interface.

The corresponding spatio-temporal diagrams, obtained by averaging both radially and azimuthally to obtain an average fluid concentration at each axial position and plotting the results as a function of time, are provided in Fig. 7 for the vertical, concentric annulus. Here, both vertical and horizontal axes in the diagrams have been normalized by the axial length $L^$⁠, which corresponds to the distance between the initial fluid interface and the outlet of the domain. In each diagram, two lines are added to indicate the approximate positions of the front and the rear of the mixing zone between the two fluids. These have been determined by tracking the axial position where the fluid concentration is 0.01 or 0.99, and are to be taken as indicative of the front and rear positions.

The main qualitative features in Fig. 6 can be recognized in these spatio-temporal diagrams. That is, at Fr = 5.1, the front of the mixing zone propagates at a near-constant speed, with a relatively smooth axial variation in the concentration field behind the front. Due to the small influence of buoyancy at this large value of imposed bulk velocity, the interface between the fluids remains well-defined. The front velocity can be estimated from the slope of the front position in the spatio-temporal diagrams, and for Fr = 5.1 a near-constant dimensionless velocity of 1.5 is found from Fig. 7. This value agrees with the peak velocity of a pressure-driven Newtonian fluid in a plane channel, whose peak velocity is 3/2 times that of the bulk velocity.

As the imposed velocity and Fr decrease, the front velocity shows some minor fluctuations due to the increased importance of buoyant stresses and the emergence of instabilities in the flow. At smaller values of Fr, more spatial and temporal variations are seen within the mixing zone as well, consistent with the radially and azimuthally averaged diagrams in Fig. 6. As observed above, the rear position of the mixing zone shows slight backflow at Fr = 0.5 and significant backflow at Fr = 0.1 due to the increased importance of buoyancy. Finally, and as expected, the largest dimensionless front velocity of the cases considered in Fig. 7 is for Fr = 0.1, i.e., the case with the lowest imposed bulk velocity. In this case, buoyancy becomes more important than in the other cases considered, resulting in both a greater, normalized front velocity and a pronounced backflow of lighter fluid.

Previous studies of density-unstable exchange flows and displacements have found displacements to exhibit macroscopic diffusion behavior under certain conditions. Alba et al. observed in pipe geometries that fully diffusive displacements corresponded to experimental conditions where the front velocity was approximately the same as the imposed bulk velocity [38]. To examine the current annular displacements, the cross-sectional average of the concentration profile at different times are plotted against a similarity variable $(x^−x^0−U^t^)/t^$⁠, where $x^0$ is the initial interface position at the start of the displacement and $U^$ denotes the imposed bulk velocity. Results for Fr = 0.1, 0.5, and 1.2 are shown in Fig. 8. The solid lines represent the solution to the 1D diffusion equation, i.e., $c=(1/2)erfc((x^−x^0−U^*t^)/2D^mt^)$⁠, where erfc() denotes the error function and $D^m$ can be viewed as a “macroscopic” diffusion coefficient. This coefficient has here been estimated using a least-squares fit to the scaled concentration profiles. Not shown in the figure is the cross-sectional average for Fr = 5.1. This case corresponds to relatively high imposed velocities and did not exhibit any tendency toward collapsing onto a macroscopic diffusion curve either. In other words, we consider none of these displacements to be particularly well represented by the diffusion-solution. This conclusion is in line with the regime classification map by Alba et al., where the diffusive regime was found for higher values of Re/Fr and/or Fr [38].

The results presented above were obtained for the special case of a perfectly vertical and perfectly concentric annulus. In this case, gravity acts exactly parallel to the axis of the annulus, and the exact centralization of the inner casing implies a uniform gap width between inner and outer walls and equal resistance to flow everywhere. Under more realistic field conditions, this perfect symmetry is broken, either due to eccentric position of the inner casing and/or due to an inclination from the perfectly vertical orientation. In this section, we will study how a minor inclination of 10 $deg$ from vertical will affect the displacement.

In Fig. 9, radial gap-averaged concentration fields are shown for Fr = 0.1, 0.5, and 1.2, while for Fr = 5.1 also azimuthally averaged diagrams are shown. Starting with Fr = 5.1, where buoyancy can be considered a perturbation to passive advection, the gap-averaged diagrams at the two points in time show that the interface now propagates faster on the lower side of the concentric annulus. That is, even at this minor inclination, the small transverse component of gravity is sufficient for accelerating dense displacing fluid toward the lower side. This trend increases in importance as Fr decreases and the imposed velocity decreases. At Fr = 1.2 and $t^=5$ s, dense fluid has been drawn toward the lower side of the annulus, similar to the case at Fr = 5.1 and $t^=6$ s, but the front has not advanced as far due to the lower imposed velocity. When comparing the concentration profiles at $t^=15$ s, shown for Fr = 0.1, 0.5, and 1.2, we observe similar gap-averaged profiles. The transverse component of gravity has been given the same amount of time to accelerate dense fluid toward the lower side of the annulus and further downward. The main difference between these three cases is that backflow is reduced as the imposed velocity increases. While the focus here is on the effect of inclination, we anticipate that, at least to a certain degree, eccentricity will have a similar effect. In a vertical, eccentric annulus, the greater gap width in the wide sector of the annulus will cause higher axial flow velocities than in the narrow sector. This effect is expected to produce a qualitatively similar trend as that observed for the mildly inclined concentric annulus,i.e., a fluid–fluid interface that elongates and propagates faster in one sector of the annulus compared to the sector at the opposite side.

In a recent study of density-unstable vertical pipe displacements, Amiri et al. developed a 1D displacement model based on the lubrication scaling [49] and the assumption that the dense fluid propagates downward through the center region of the pipe, with the lighter fluid occupying the annular region between the center core and the pipe wall [39]. The 1D lubrication model represents a balance of viscous and buoyant stresses, and can therefore be useful to make predictions concerning onset of backflow or front velocities for non-inertial displacements. In the simulation results presented in Fig. 6, the azimuthal averages (2) suggest a similar tendency for propagation of dense fluid through the center of the annular gap. This tendency is most visible at the highest values of Fr in the figure. This suggests that a similar 1D displacement model could be applicable for vertical, concentric annuli as well, in the limit where flow is primarily in the axial and radial directions. Below, we use the 1D lubrication model for a vertical concentric annulus, treated as a plane vertical channel, that is listed in the appendix of Ref. [39]. Some details of the derivation are included below and in the appendix for completeness.

To see the effect of χ on the overall displacements within the lubrication model, Fig. 10 shows the time evolution of the interface between the two fluids at different values of χ. At relatively small values of χ, the viscous stress dominates over buoyant stresses, and there is a net downward flow of both fluids. A shock develops along the center of the channel (corresponding to low values of c in Fig. 10) for positive values of χ. At larger values of χ, the profile steepens into a shock toward the wall (at higher values of c), and at the critical value χ* = 117.8, the shocks at the wall become stationary. This critical value of χ* is in agreement with the value listed by Amiri et al. [39], as is the extent of the stationary layer in Fig. 10. At even larger values of χ, there is a net backward (upward) flow of lighter fluid along the walls, as seen for the case χ = 150 in the figure.

The solutions to these equations are shown in Fig. 11, along with the front velocities from 3D CFD simulations which have been determined from the spatio-temporal diagrams. The error bars attached to the data points from the vertical CFD simulations represent the variation in computed front velocity when the volume fraction is varied from 0.01 to 0.05. Clearly, the largest variation is associated with the most buoyant simulation, which is to be expected as these are the more unstable cases considered. Error bars have not been attached to the inclined simulations, as the choice of volume fraction for determining the front position did not affect the numerical values in Fig. 11 for the inclined annulus. The front velocity predicted by the simplified model is in reasonably good agreement with the CFD simulations in the vertical and concentric annulus, indicated by the filled points. Excellent agreement is obtained at the lowest values of χ which corresponds to the simulations with Fr = 1.2 and Fr = 5.1. These are the more stable displacements, where the averaged fluid concentration profiles agree with the simplified model assumptions of a well-defined, single-valued interface and where viscous stresses are significant compared to the effect of buoyancy. The 1D model under-predicts the front velocity at greater values of χ, but the agreement is still considered reasonable. Also indicated in the figure are the front velocities from 3D CFD simulations of the inclined annulus. For the three lowest values of Fr, these are indicated by empty circles in Fig. 11 and the dashed line is simply a guide to the eye, corresponding to a linear least-squares fit to the data. As expected, by breaking the symmetry of the original model, in this case by inclining the annulus, the front speed increases significantly compared to the symmetric base case. A similar increase in front velocity due to inclination has also been observed both for exchange flows and displacement flows in pipe geometries [33,35,38].

Concerning the front velocity prediction by the 1D model, we note that this is a function that is nearly linear in χ. This is due to the width of the shock, c*, which, for values of χ smaller than approximately 50, increases rapidly toward a magnitude close to 0.25. This is indicated by the dashed-dotted curve in Fig. 11. As χ grows, c* slowly approaches a value of approximately 0.265. Thus, since c* is either small (for small values of χ) or close to a constant value of 0.25–0.26 for larger values of χ, Eq. (8) will be close to linear in χ.

Finally, concerning backflow of lighter fluid, the 1D model predicts the transition between instantaneous displacement and backflow to occur at χ* = 117.8. For this critical value of the dimensionless parameter, the shock adjacent to the wall (where the light fluid is assumed to be located) is static. As shown in Fig. 10, higher values of χ would produce backflow according to this model. Setting Re/Fr = 55, as per the simulations studied here, this corresponds to a critical value of Fr ≈ 0.94, below which the 1D model predicts backflow. This is in agreement with the spatio-temporal diagrams for the vertical, concentric annulus in Fig. 7, where Fr = 0.1 and Fr = 0.5 exhibit backflow, while Fr = 1.2 and Fr = 5.1 do not. While this simplified model provides reasonable agreement with the 3D CFD simulations for the vertical, concentric annulus, it is clear from the above that effects of inclination (or eccentricity) would not be correctly reflected in the 1D model that neglects azimuthal variations. A gap-averaged 2D Hele–Shaw model, such as that of Bittleston et al. [51] and Maleki and Frigaard [52], is expected to be a significantly better approximation in these cases and could be the topic of future studies.

The field case reverse circulation operation that served as the starting point for defining this numerical study exhibited a well-defined top-of-cement relatively few intrusions of foreign material [26]. These observations, along with the relatively low permeability of the two retrieved sandwich joints, suggest that the reverse circulation cementing operation was largely successful, at least locally from the regions where the joints were retrieved. The 3D numerical simulations of this case presented above (Fr = 5.1) suggest that the destabilizing effect of the density difference is dominated by viscous stresses at this imposed velocity. In the perfectly symmetric vertical and concentric case, this is evidenced by the sharp interface between the fluids that appeared to be advected with the flow, and that the interface propagation speed is approximately the same as the peak velocity in the center of the channel. The main effect of the density difference at Fr = 5.1 in Fig. 6 is probably the flattening of the front of the interface, seen in the right-most diagram in Fig. 6, similar to the center-channel shock predicted by the 1D lubrication model. The effect of density difference at Fr = 5.1 is more visible in the inclined annulus, where the transverse component of the gravitational acceleration promotes flow of dense fluid toward the lower side of the annulus. By decreasing the imposed bulk velocity, viscous stresses decrease and the displacements turn increasingly unstable. The result is increased tendency for backflow and a front velocity that propagates significantly faster than suggested by the imposed flow alone. At field conditions, Fr = 5.1, viscous stresses dominate buoyancy and result in a net downward flow with no apparent backflow of the lighter fluid. Thus, the field case appears to be successful since viscous stresses are allowed to dominate the destabilizing buoyant stresses, through a combination of small density difference, high imposed velocities, and high fluid viscosity. Based on the above observations and the regime classification suggested by Alba et al., we categorize Fr = 1.2 and Fr = 5.1 as viscous and instantaneous displacements at $Re/Fr=55$ and At = 0.022. The simulations with Fr = 0.1 and 0.5, corresponding to greater values of χ, are categorized as non-instantaneous (i.e., with backflow) and inertial due to their more significant azimuthal instabilities and poorer agreement with the front velocity predicted by the lubrication model.

This study has taken a specific reverse circulation operation as starting point for defining geometry and fluid parameters. Several simplifications have been made, such as setting the viscosity of both fluids to the same, constant value, and treating the annulus as perfectly concentric. Forces acting on the casing will in practice make the casing off-centered or eccentric, even in near-vertical hole sections [53]. This was also the case for the retrieved sandwich joints, shown in Fig. 2. The local geometry will affect both the displacement and frictional pressure losses [54]. Furthermore, the study has considered simple, constant viscosity fluids, and does therefore not address yield stress and shear thinning effects on density-unstable displacements [40]. While the Newtonian character of the cement slurry was based on actual measurements of what is perceived to a comparable slurry formulation, the displaced fluid (i.e., the drilling fluid and/or spacer fluid) will normally be characterized by a nonlinear viscosity function and exhibit a non-zero yield stress. In a highly eccentric annulus, such fluids will become harder to mobilize from the narrow side [15] and this may exacerbate the situation by further promoting channeling of dense fluid along the wide sector of the annulus. Clearly, more studies are required to explore different viscosity hierarchies and more realistic, eccentric annuli. The combination of eccentricity and inclination may in particular result in interesting reverse circulation displacements. Assuming that gravity causes the inner casing to be shifted toward the lower side in an inclined wellbore, the eccentricity should now tend to offset the effect of inclination, since the gap width in the lower side is narrower than at the top side of the annulus. Future work in this direction should attempt to study these more realistic conditions, as well as explore different displacement categories, similar to Ref. [38] for pipe geometries.

Primary cementing of casings is normally performed by injecting spacer fluids, washes, and cement slurry down the well inside the casing. The drilling fluid occupying the annular space behind the casing will therefore be displaced from bottom and up toward the surface, and normally with a density-stable hierarchy of fluids. Design rules based on experience and experimental results have been developed to facilitate the job design, as well as fundamental insights based on analyses of the displacement mechanics. While conventional circulation is by far the more common approach to cementing, under certain conditions, reverse circulation cementing may be preferred or required. Reverse circulation cementing can reduce the circulating pressures at the bottom of the hole, enabling higher displacement rates and denser cement slurries. However, this approach also reverses the order of the fluids in the annulus, placing the denser cement slurry above lighter fluids to be displaced, resulting in a mechanically unstable configuration and potentially more mixing between fluids and a larger risk of slurry contamination.

This study has used an actual reverse circulation cementing operation as starting point for defining fluid properties and injection velocity. Results from fully 3D CFD simulations suggest that the relatively high imposed velocity that was probably used during the operation generated sufficiently high viscous stresses to largely suppress the destabilizing effects of the density difference. This agrees at least qualitatively with the conditions of the cement, which have previously been characterized both visually, by surface re-logging and permeability measurements. This initial study has focused on the most symmetric case of a vertical, concentric annulus, and an annulus with a mild inclination from vertical. It is anticipated that both increased inclination as well as eccentricity can result in more adverse conditions for reverse circulation cementing, such as tendencies for stratification or channeling of the dense fluid and a longer mixing zone. At lower imposed bulk velocities, buoyancy becomes increasingly important, resulting in instabilities near the interface and an increased tendency for backflow of lighter fluid. Results of this initial study suggest that a higher injection rate can act to suppress this instability, at least in close-to vertical annuli that are well centered. A classification of the displacement simulations carried out in this study suggests instantaneous (i.e., no backflow) and viscous displacements for the highest values of Fr considered, i.e., Fr = 1.2 and 5.1. At lower values of Fr, the displacements approach the exchange flow limit, and here the displacements are non-instantaneous and exhibiting more pronounced azimuthal and radial instabilities. These are categorized as inertial, in line with the regime classification presented by Alba et al. [38]. Future work will address more realistic eccentricities and fluid viscosities and attempt to categorize displacements according to the governing dimensionless numbers.

This work was supported by the Research Council of Norway through project number 294815 and performed utilizing computational infrastructure provided by UNINETT Sigma2—the National Infrastructure for High Performance Computing and Data Storage in Norway.

There are no conflicts of interest.

The data and information that support the findings of this article are freely available. The authors attest that all data for this study are included in the paper. Data provided by a third party are listed in Acknowledgment.

c =

volume fraction of the displacing fluid

q =

flux function

$g^$ =

gravitational acceleration, L T⁻², m/s²

$h^$ =

distance from center of annular gap to fluid interface, L, m

$p^$ =

pressure, M L⁻¹ T⁻², Pa

$t^$ =

time, T, s

$v^$ =

velocity, L T⁻¹, m/s

$x^$ =

axial position, L, m

$L^$ =

axial length, L, m

$D^mol$ =

diffusion coefficient, L² T⁻¹, m²/s

$D^m$ =

“macroscopic” diffusion coefficient, L² T⁻¹, m²/s

$T^*$ =

characteristic time scale, T, s

$D^h*$ =

characteristic length scale, L, m

At =

Atwood number

$D^i,D^o$ =

inner (i) and outer (o) annular diameter, L, m

Fr =

densimetric Froude number

Pe =

Péclet number

Re =

Reynolds number

u, v =

dimensionless velocity

$μ^$ =

viscosity, M L⁻¹ T⁻¹, Pa s

$ρ^$ =

mass density, M L⁻³, kg/m³

ϕ_i =

integer fluid identifier (±1)

χ =

ratio of buoyant and viscous stresses

In this section, the 1D displacement model based on a lubrication scaling will be derived. The derivation assumes viscous displacements, where the dense displacing fluid occupies the center of the vertical, concentric annulus, with the lighter adjacent to the walls of the annulus. A similar model was proposed by Flumerfelt [50] for vertical, density-stable displacements. The assumed fluid configuration is illustrated in Fig. 12, where it is assumed that the fluids are separated by an interface at distance $y^=h^(x^,t^)$⁠, from the center of the annulus. It is assumed that $h^$ is single-valued and symmetric about the center of the annular gap.

The governing equations are simplified by invoking a lubrication scaling that assumes that axial variations of the streamlines occur over a length scale $L^*$ much larger than that of the radial gap, $D^h*$ [49]. A similar model has been developed for miscible displacements in a vertical pipe by Amiri et al. [39].