Reverse Circulation Displacement of Miscible Fluids for Primary Cementing

The construction of wells for petroleum production, geological carbon storage, or geothermal energy recovery proceeds in stages by first drilling the wellbore to a certain depth, and next running and cementing a casing string inside the hole. The operation of cementing the casing is referred to as primary cementing, and involves displacing the drilling fluid from the annular space outside the newly run casing, and replacing it with a cement slurry. Once placed, the cement slurry should harden to form an annular barrier to protect the casing string from mechanical loads and corrosive formation fluids, and provide zonal isolation along the wellbore [1]. After cementing the casing, the well can be extended by drilling out the cement at the bottom of the casing, pressure testing of formation, and annular cement [2], followed by drilling of new formation.

Displacement of drilling fluid from the annular space behind the casing is normally performed by pumping cementing fluids such as washes, spacer fluids, and cement slurry down the well inside the casing to be cemented. This cementing design is referred to as conventional circulation cementing. Wiper plugs are typically used to separate fluids while flowing down inside the casing. Once the cementing fluids exit the casing at the bottom of the well, and flow up toward the surface in the annular space behind the casing, the washes and spacer fluids should ensure minimal mixing between drilling fluid and the cement slurry [1]. Recommended practices for effective fluid displacement typically involve centralization of the casing and stabilization of the fluid displacement by maintaining density and viscosity hierarchy between successive fluids [3].

An alternative primary cementing strategy is to inject the cementing fluids directly to the annulus that is to be cemented [4–6]. This design is often referred to as reverse circulation cementing which can be performed when the infrastructure allows direct access to the annulus from surface [7], such as in most onshore wells and offshore platform wells. Reverse cementing was motivated initially by the lower circulating pressures toward the bottom of the hole compared to conventional cementing [8], allowing for denser cement slurries across weak and depleted formations, higher flowrates during displacing drilling fluid. Further, the use of retarders can be reduced compared to conventional cementing, since the slurry is now placed directly into the annulus [8–10]. This saves time since the cement will cure faster when less retarders are added and results in improved compressive strength development of cement particularly in the upper hole sections [9,11]. 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 [12].

As most cement slurries used in primary cementing are denser than typical drilling fluids, reverse circulation cementing results in a density-unstable fluid configuration in the annulus during mud displacement and placement of cement slurry. This is linked to the Rayleigh–Taylor instability [13], and can result in increased intermixing between the fluids as well as periods of free-fall unless cementing fluids are pumped sufficiently fast [14]. The recent studies of Eslami et al. [5] and Skadsem and Kragset [6] have used respectively experimental and computational methods for studying buoyant reverse circulation displacements. In the experimental study of Eslami et al., effects of density difference, imposed flow, centralization, fluid viscosity, and reciprocation of the inner string were all considered. The displacements were examined primarily by imaging and having the displaced fluid colored by black ink. By linking the recorded light intensities to mixtures of the displaced and displacing fluids [15], the images were mapped to the spatiotemporal fluid concentration field, which in turn provides information concerning front velocities and the structure of mixed fluid regions. Their results show that buoyant displacements become less effective and more unstable as the displacing fluid is made less viscous [5]. Eccentricity was found to be significant in the case of viscoplastic fluid displacement, but less so for Newtonian fluid configurations.

The computational study by Skadsem and Kragset focused on buoyant, reverse circulation displacements in a concentric and vertical annulus [6]. The study was motivated by pressure and fluid leakage measurements conducted on retrieved sections from an abandoned North Sea well that had been cemented using the reverse circulation method in the mid-1980s [4,16,17]. The computational study was limited to a fixed fluid hierarchy, i.e., iso-viscous and Newtonian fluids, with a fixed density difference, and varying imposed flow velocities. Using fluid properties and flowrate based on the original job report, the computational study showed that the largest imposed flowrate resulted in a relatively steady displacement where large viscous stresses largely suppressed destabilizing buoyant effects. At lower flow velocities, more mixing and a tendency toward backflow of light fluid was observed. While the main thrust of the study was the results obtained using fully three-dimensional computational simulations, it was also shown that a one-dimensional displacement model based on lubrication scaling could predict front velocity reasonably well for the perfectly vertical annulus. Even a modest inclination from vertical (10 deg) was found to significantly increase the front velocity. Inclination results in a preferred orientation of the fluids (dense fluid toward the low side, light fluid toward the high side), which results in increasing front velocity [6].

Buoyant exchange flows and buoyant displacements have been studied in the past using other geometries, such as inclined pipes and ducts. The exchange flow of two miscible fluids were studied by Debacq et al. in a vertical, transparent pipe [18,19]. By coloring one of the fluids, and by calibrating video images to equivalent fluid concentration, cross-sectional average fluid concentration profiles were generated for each experiment. The experiments showed that the mixing region between the two fluids evolved according to a diffusive transport mechanism, characterized by a macroscopic diffusion coefficient that was estimated to be orders of magnitude larger than typical molecular diffusion coefficients [18]. Identification of diffusive mixing was made by scaling the axial position, $z^$⁠, of the concentration profiles according to the diffusive similarity variable, $z^/t^$⁠, where $t^$ denotes the time since start of the experiment. For diffusive displacements, axial concentration profiles scaled in this manner should approach a single curve, i.e., the similarity solution to a one-dimensional diffusion equation, characterized by a diffusion coefficient [18]. The buoyant mixing of the two fluids was found to be diffusive over a large range of viscosity ratios and tube diameters, and could be characterized as either convective or turbulent diffusive [19]. A third, stable counter-flow regime was also observed, predominantly for exceedingly small density differences and in small diameter tubes [19]. Imposing a fixed velocity to the fluids result in buoyant displacement flows, which have been studied by e.g., Taghavi et al. [20,21] and Alba et al. [22]. An extensive catalog of experiments have demonstrated different displacement regimes, involving net upward flow of the lighter fluid (backflow), inertial, viscous, or macroscopically diffusive displacement flows [22], depending on pipe inclination and primarily two dimensionless numbers, i.e., the Reynolds number and the densimetric Froude number [22].

We present a combined experimental and computational study of reverse circulation displacement of the annulus, where the range of displacement conditions is extended beyond our previous computational study [6]. A transparent laboratory-scale annulus test rig has been designed and constructed for this purpose, and a mirror arrangement is used to image the fluid interface simultaneously from three different directions. We focus on iso-viscous displacements in a vertical and concentric annulus, and vary the density hierarchy among the fluids to study both stable and density-unstable displacement conditions. Numerical simulations have been performed to further explore the buoyant displacements using the openfoam computational platform and a numerical scheme previously used to simulate buoyant exchange flows [23].

Displacement experiments were conducted using a dedicated annular flow loop consisting of an outer transparent pipe with inner diameter 40 mm and an inner opaque pipe with outer diameter 25 mm. Primary cementing of casings for well construction involves displacements in relatively narrow annular spaces between successive casing strings [24]. Commonly used combinations of casing strings and hole sizes result in ratios of inner to outer annulus diameters (annular aspect ratios) corresponding to $κ=D^inner/D^outer∼0.7−0.8$⁠. The corresponding aspect ratio for the annular test rig is 0.63, which is somewhat less than that of common wellbore geometries. The slightly wider annular geometry in the experimental setup enables improved imaging of the displacement flows, while still retaining the relatively narrow radial clearance between the two pipes. In all experiments, the inner pipe was centralized inside the outer pipe, forming a concentric annulus, and the annular test section was aligned vertically, so that gravity acted parallel to the imposed flow direction, which was from the top. The position of the inner pipe was controlled by fitting specially designed rigid centralizers, produced by additive manufacturing, to the top and bottom of the inner pipe. A manual ball valve was used to separate the fluids before the start of each experiment. The annular test section is 1500 mm long, and extends from the position of the ball valve (corresponding to the initial interface position) to the bottom outlet of the flow loop. The length of the test section is thus 100 hydraulic diameters. When presenting the experimental results below, the initial interface position is taken as the origin in the axial direction, $z^=0$⁠. We note that the design of the flow loop prevented image capturing above the initial interface position. This is not a limitation for stable reverse circulation displacements, where the displacing fluid is moving downward, parallel to the imposed flow direction. In density-unstable experiments, however, the rig design implies that we were unable to visualize backflow of light, displaced fluid above the initial interface position. In our numerical simulations, we do capture the flow field above the initial interface position.

Black dye was added to the displacing fluid to enable tracking of the fluid interface and the mixing region during displacement experiments. A Nikon D5500 camera operating at 60 frames per second was used to capture images of the experiments. The camera was placed in front of the test rig, and a mirror arrangement allowed simultaneous imaging of the annulus from three different directions, as shown in Fig. 1. A rectangular tank filled with glycerol was attached outside the annular test section to improve imaging of the fluids. An Endress+Hauser Promag 53 flowmeter was fitted to the outlet of the test section to control the imposed flowrate in displacement experiments. Finally, an absolute pressure sensor was installed at the bottom of the test section. All experiments reported in this study were gravity-driven. At the inlet to the annular test section, an arrangement consisting of an overflow tank ensured a constant vertical liquid level throughout each experiment. A manual choke valve at the outlet was regulated to produce different magnitudes of imposed flow, or kept closed in the case of buoyant exchange flow experiments.

Tap water and brine were used as the miscible test fluids in this study. Sodium chloride was added and dissolved in tap water to produce the brine, which was the denser of the two fluids. Fluid densities were measured using a DMA 4500 M density meter. The experimental program encompasses density-stable, iso-dense, and density-unstable displacements, and three different imposed velocities, as outlined in Table 1. Here, and in what follows, $ρ^1$ refers to the density of the displaced fluid, and $ρ^2$ to the density of the displacing fluid. The ratio of the density difference to the average density is expressed by the Atwood number, $At=(ρ^1−ρ^2)/(ρ^1+ρ^2)$⁠. Finally in Table 1, $U^0$ denotes the magnitude of the imposed (bulk) flow velocity.

No viscosifier was used in any of the fluids, so the fluid pairs in these experiments are considered iso-viscous. As specified in the table, exchange flow experiments with no imposed flowrate were also performed for the buoyant cases (series 3 and 4). Experiment series 1 and 2 in Table 1 correspond to iso-viscous density-stable, and iso-viscous density neutral displacements, respectively. Reverse, or downward, density-stable displacements have been considered experimentally and theoretically by Lajeunesse et al. in a rectangular channel or Hele–Shaw geometry [25]. For iso-viscous fluids and the Atwood number and imposed flowrates associated with series 1 listed in Table 1, their theoretical study showed that the interface will take on a steady shape that exhibits an internal shock (jump discontinuity), and with a foot advancing ahead of the shock through the center of the gap [25] (i.e., domain 2 in Ref. [25]). The iso-dense and iso-viscous displacements, corresponding to series 2, are equivalent to one and the same fluid displacing itself. As such, the fluid-fluid interface is expected to advect axially according to the concentric annulus velocity profile, and with the front of the interface advancing at a speed slightly greater than $1.5U^0$⁠. Finally, series 3 and 4 represent reverse circulation density-unstable exchange flows and displacements with imposed velocity. These represent cases of the classical Rayleigh–Taylor instability where gravity will act to destabilize the interface.

Relevant dimensionless numbers can be identified by scaling the governing equations as per Ref. [6]. To this end, the imposed bulk velocity $U^0$ is taken as the characteristic velocity scale and the gap width of the annulus $D^h*=(D^o−D^i)/2$ as a characteristic length scale. Using the viscous stress $μ^U^0/D^h*$ as relevant scale for the pressure, we identify the Reynolds number, $Re=ρ¯^U^0D^h*/μ^*$⁠, with $ρ¯^$ being the average mass density of the fluids, the densimetric Froude number $Fr=U^0/Atg^D^h*$⁠, and the Péclet number $Pe=U^0D^h*/D^mol$ as the main dimensionless numbers. As per our previous study, we assume Pe → ∞ so microscopic diffusion can be neglected in Eq. (3).

For the numerical solution of the governing equations, the volume-of-fluid solver twoLiquidMixingFoam is used. This solver is specialized for two-phase, incompressible, and miscible fluid flows. Second order discretization schemes are used for the spatial terms and time discretizations. For the boundaries, 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. To stabilize the solution an adjustable time-step with a maximum Courant number of 0.5 is used. The computational grid is a structured grid, stretched toward the walls. It consists of 120 grid cells in the azimuthal direction, 22 cells in the radial direction, and 900 cells in the axial direction. The computational grid is shown in Fig. 2.

In the following sections, results from the experiments listed in Table 1 will be presented mainly in the form of spatiotemporal diagrams. These diagrams show the cross-sectional average concentration, $c¯$⁠, as a function of position and time, and are practical for comparing the averaged concentration profiles across different experiments, and gauge the speed of the advancing front of the interface as well as the overall stability of the displacement.

Frontal view images of density-stable displacement (series 1) at an imposed velocity of 12 mm/s are shown in Fig. 3(a). Consecutive images are taken approximately 13 s apart, and show a symmetric and steadily advancing interface. Shown in the right panel of Fig. 3 is the corresponding spatiotemporal diagram of this experiment. The displacing fluid corresponds to $c¯=1$⁠, and displaced fluid is taken as $c¯=0$⁠. The superimposed dashed line is indicative of the interface position, and suggests that the interface advances downstream at a fixed speed, approximately given by the imposed bulk velocity of 12 mm/s. The sudden transition from $c¯=0$ to $c¯≈1$ suggests a minimal mixing region between the fluids and a practically horizontal interface that is translated due to the imposed velocity. These observations are qualitatively similar to comparable results presented by Lajeunesse et al. [25]. Similar results are also found for the other density-stable displacements (series 1), but with more rapidly advancing interface due to the greater imposed velocities (25 mm/s and 33 mm/s).

Experimental spatiotemporal diagrams for the iso-dense displacements (series 2) are shown in Fig. 4 for three different values of the imposed flow velocity (corresponding to series 2 in Table 1). Compared to the density-stable case considered in Fig. 3(b), the iso-dense displacement in Fig. 4(a) exhibit a mixing zone that spreads as the fluids are transported axially, as expected. Further, comparing Figs. 3(b) and 4(a), the front is seen to advance more rapidly in the latter experiment, and approximately 1.5 times faster than the density-stable experiment. This is in accordance with an interface that advances according to the laminar velocity profile in a concentric annulus, where the peak value of the velocity is approximately $1.5U^0$⁠.

As outlined in Table 1, two different density ratios have been studied for density-unstable experiments, corresponding to At = 0.003 and At = 0.02. Exchange flow experiments, with no imposed flow, and displacement experiments with imposed flow have been considered.

Spatiotemporal diagrams for the two exchange flow experiments are shown in Fig. 5. The flow is here driven by buoyancy alone, similar to the classic Rayleigh–Taylor instability, but now studied within the confined space of the vertical, concentric annulus. Compared to the stable and the iso-dense configurations shown above, mixing and fluctuating average concentration profiles are seen in the spatiotemporal diagrams for these density-unstable displacement experiments. A relatively well-defined front, indicated by the dotted curve in Fig. 5, is seen to advance at an approximately constant speed (corresponding to the slope of the curve). Fluctuating concentration levels follow in the wake of the interfacial front position. The slope of the curve that tracks the front of the mixing zone indicates the front velocity of the displacement. At late times, we observe a relatively constant front velocity for At = 0.003, in Fig. 5(a). The larger density difference, represented here by At = 0.02, results in a more unsteady advancing front velocity. The dashed curve in Fig. 5 is a straight line fitted to the late-time front position, and its slope is taken as the approximate front velocity, $U^f$⁠, for the two cases. We find $U^f≈15mm/s$ for At = 0.003, and $U^f≈29mm/s$ for At = 0.02. The larger density difference causes the faster front propagation, as expected. A natural velocity scale for buoyant flow is $U^g=Atg^D^h$⁠, which we evaluate to 21 mm/s and 54 mm/s for At = 0.003 and At = 0.02, respectively. Thus, the ratio $U^f/U^g$ is approximately 0.7 for At = 0.003 and 0.5 for At = 0.02.

We next turn to buoyant displacements with an imposed flow, focusing first on series 3 in Table 1. In Fig. 6, the left panel shows front view images captured at 4 s intervals, while the right panel is the corresponding spatiotemporal diagram showing the evolution of the cross-sectional average concentration. This displacement involves the downward propagation of dense fluid, combined with the tendency for the light, displaced fluid to flow upward, in the opposite direction of the imposed flow. As pointed out above, the front velocity in the case of buoyant exchange flow (no imposed flow velocity) was measured to approximately 15 mm/s for At = 0.003. In Fig. 6, the imposed velocity is 12 mm/s and the front velocity is approximately $U^f≈28mm/s$⁠. That is, the front velocity can in this case be approximated by the sum of the exchange flow front velocity and the imposed flow velocity.

In Fig. 7, a comparison is provided of spatiotemporal diagrams for three different imposed velocities and At = 0.02. The time scale of the experiments and the interface front position suggest that the interface advances with an increasing speed as the imposed flow velocity increases, as expected. One may also observe a slight stabilization of the displacement as the imposed flow velocity increases in the form of suppressed fluctuations within the mixing zone, and a more complete displacement in the rear. As discussed in connection with Fig. 5, the speed at which the front of the mixing zone propagates downward may be inferred from the spatiotemporal diagrams by tracking the (approximate) front position in time. For the buoyant exchange flows with no imposed flow, the ratio of front velocity to characteristic buoyant velocity was found to be 0.5–0.7. Considering the combined effects of an imposed flow, $U^0$⁠, and buoyancy, we estimate the front speed from the late-time front position evolution in the spatiotemporal diagrams, and consider the scaling $U^f∼U^0+0.5U^g$ in Fig. 8. For the largest density difference, At = 0.02, this scaling produces a reasonably accurate estimate of the front speed for the exchange flow experiment and the three imposed velocities considered. The front velocity is slightly larger than $U^0+U^g/2$ for the lower density difference.

As shown previously for miscible buoyant exchange flows and displacements in inclined tubes [18,22,26], the cross-sectional average concentration profile can take on a self-similar shape, corresponding to an axial diffusive spreading of the averaged concentration profile. That is, the center of the mixing region (⁠$c¯≈0.5$⁠) translates downstream with the bulk imposed velocity, and the extent of the mixing region scales according to a diffusion equation, with an effective “macroscopic” diffusion coefficient. To investigate whether the current experiments exhibited mixing regions that spread according to a self-similar shape, late-time concentration profiles have been extracted from each experiment, and plotted as function of the scaled axial position, $(z^−U^0t^)/t$⁠, as per Alba et al. [22]. As suggested above, this choice of scaled coordinate is motivated by the self-similar solution to the one-dimensional diffusion equation, which is in the form of a complementary error function of the variable $z^/4D^t^$⁠, with $D^$ the diffusion coefficient. To account for advective transport due to the imposed flow velocity, the axial coordinate is shifted by $U^*t^$⁠, as indicated above. Experiments where the late-time concentration profiles approximately collapse to the self-similar solution of the diffusion equation, may reasonably be categorized as diffusive, with a mixing zone that increases in length according to a “macroscopic” diffusion coefficient, $D^$ [22].

In Fig. 9, the scaled late-time concentration profiles are considered for experiments with At = 0.02, corresponding to the largest density difference considered. As explained, the axial coordinate is scaled according to self-similar diffusive behavior, and the superimposed solid curve corresponds in all cases to the similarity solution, i.e., $c¯=(1/2)erfc[(z^−U^0t^)/4D^t]$⁠, with erfc() denoting the complementary error function [27]. The value of the ”macroscopic” diffusion coefficient, $D^$⁠, is determined manually for each case. Arguably the best fit to the similarity solution is found for $U^0=25mm/s$⁠, shown in Fig. 9(b), where the scaled experimental curves intersect $c¯=0.5$ at $x^=0$⁠, as required by the similarity solution. The correspondence between the two other cases and a similarity solution is less clear and will be further explored when assessing the numerical simulations below. Results for the smaller density difference, corresponding to At = 0.003 and series 3 in Table 1, are shown in Fig. 10. Fair agreement is found for $U^0=12mm/s$ and $U^0=25mm/s$⁠, Figs. 10(a) and 10(b), while the case $U^0=33mm/s$ deviates from the similarity solution, in particular close to the rear (upper) part of the mixing region.

To gain further insight into the displacements studied experimentally and discussed above, and in particular to assess the degree of backflow and possible self-similar spreading of the mixing region, CFD simulations of the density-unstable cases in Table 1 have been performed. In Fig. 11, simulation results are presented for the same case as shown in Fig. 6 above. While the experimental setup does not permit visualization of backflow, the numerical simulations do. As shown in the spatiotemporal diagram in Fig. 11(b), the initial interface position is located 0.5 m away from the inlet point in the model geometry, and 1 m above the outlet point. Included in the spatiotemporal diagram in Fig. 11(b) as the dotted curve, is the experimentally determined front position (also shown in Fig. 6(b)). The curve is included for comparison purposes, and we observe good quantitative agreement between experimental and numerical predictions for the front position. The spatiotemporal diagrams and the “front view” presentation of the displacement flow is also qualitatively similar, as can be seen when comparing Fig. 11 and Fig. 6.

In Fig. 12, spatiotemporal diagrams based on 3D CFD simulations are shown for At = 0.02, corresponding to the experimental results in Fig. 7. Also here, the front position of the mixing region, as determined in the corresponding experiments, are shown as the dotted curves for comparison with the numerical results. First of all, fair quantitative correspondence between front positions are seen, particularly for the lowest imposed bulk velocity, in Fig. 12(a). As above, features within the mixing zone are qualitatively similar to the experimental results in Fig. 7; the front velocity increases with increasing imposed velocity, and increasing the imposed velocity tends to impart a stabilizing effect in that backflow is suppressed. For At = 0.02, all displacements exhibit some backflow, i.e., presence of lighter, displaced fluid above the initial interface position over the course of the simulations.

Now, Fig. 13 illustrates late-time concentration profiles from CFD simulations with the axial position scaled according to a diffusively spreading interface, $(z^−U^0t^)/t^$⁠. The superposed solid curve is the error function solution to a one-dimensional diffusion equation with a manually fitted diffusion coefficient. For comparison, the dashed curve corresponds to the diffusive scaling of the experiments, from Fig. 9. Interestingly, the numerical simulations indicate relatively good correspondence between the scaled interface profiles and the similarity solution, more so than the experimental results in Fig. 9. When comparing the fitted similarity solutions, shown as the solid and dashed curves, the agreement is considered fairly good. The fitted macroscopic diffusion coefficients are shown in Fig. 14 for the displacements corresponding to At = 0.02. We observe that simulations consistently predict larger macroscopic diffusion coefficients, indicative of more rapid growth of the mixing region during displacement than the experiments.

In their study of buoyant displacement flows in tilted and in vertical tubes, Alba et al. found diffusion coefficients of the order of 10⁻³ m²/s, and that the diffusion coefficient increased with (i) increasing imposed velocity, and (ii) increasing inclination angle from the vertical [22]. Considering the values for the diffusion coefficient in Fig. 14, we find numerical values to be of the same order of magnitude as that found for vertical tubes [22]. As pointed out in previous studies, the numerical values of the diffusion coefficient are orders of magnitude larger than that of molecular diffusion, and this suggests that the main transport mechanisms are advection and axial dispersion (and not molecular diffusion between the fluids). The current experimental and numerical results suggest an increasing trend in the macroscopic diffusion coefficient with increasing imposed velocities, as per Ref. [22], but this is based on the relatively few observations shown in Fig. 14. The diffusion coefficients based on the simulations are not conclusive when it comes to this trend, so additional simulations and/or experiments should be performed to clarify whether trends exist also for the vertical and concentric annulus.

Finally, possible macroscopic diffusion behavior in the numerical simulations corresponding to At = 0.003 are considered, as shown in Fig. 15. Interestingly, the collapse onto a diffusive similarity solution is now not achieved for all values of $c¯$⁠, and this is particularly evident for the two largest imposed velocities, shown in Figs. 15(b) and 15(c). Similar non-diffusive behavior has been observed by Alba et al. [22], who went on to characterize the non-diffusive displacements as either inertial or viscous, depending on whether the displacement flow displays evident instabilities or not. Since all density-unstable exchange and displacement flows in this study (series 3 and 4 in Table 1) exhibit instabilities, we choose to classify the displacements with At = 0.003 as non-diffusive and inertial based on the numerical results. Similarly, the displacements with At = 0.02 shown in Fig. 13 are categorized as diffusive, again primarily based on the results of the numerical simulations.

This classification can be summarized in Fig. 16, using the same notations and symbols as in Alba et al. [22]; that is, points with superposed squares (circles) represent inertial (viscous) displacements while the unadorned points are diffusive displacements. Further, the displacements corresponding to the larger Froude numbers (Fr > 1) are instantaneous, i.e., there is no discernible backflow during the simulations. In Fig. 16, the diffusive displacements at Re/Fr = 821 correspond to At = 0.02, while the points at Re/Fr = 341 correspond to At = 0.003. Also shown are observations corresponding to Re/Fr = 55, which are taken from the previous numerical study of Skadsem and Kragset [6], and shown with the classification suggested in that study. The two straight lines in Fig. 16 correspond to χ = 2Re/Fr² = 117.8, which marks the criterion for onset of backflow in the viscous limit (see Ref. [6]), and Re/Fr = 500 − 50 Fr, which Alba et al. used to discern diffusive from inertial displacements, [22]. As such, the classification of the current vertical and concentric annular displacements are considered to be in accordance with the displacement regimes developed for buoyant displacements in pipes [22]. Although the current data set for the vertical and concentric annulus agrees well with the previous regime map developed for pipe displacements, additional studies should seek to confirm this view, and to explore effects of eccentricity and inclination on the displacement regime classification for buoyant displacements in annuli.

In summary, we have conducted a combined experimental and computational study of reverse circulation displacements in a vertical and concentric annulus. We have focused on the effects of density difference between the two fluids, covering density-stable, neutral, and unstable conditions, and potential stabilizing effect of the imposed flow velocity. We limit the study to iso-viscous fluid pairs. The density stable and density neutral experiments are in qualitative agreement with expectations, as the fluid interface is symmetric and advected steadily downstream. The density neutral fluid pair corresponds essentially to a single fluid displacing itself, and we find good agreement between the propagation speed of the mixing region and the peak magnitude of the concentric annulus velocity profile.

More relevant for cementing flows and reverse circulation displacements are the density-unstable conditions, corresponding to the displacement of lighter drilling fluids by denser cementing fluids. Density-unstable conditions result in more rapid propagation speed of the interface front, and a growing mixing region between the fluids with fluctuating concentration levels. Although not directly observed in experiments, our numerical simulations have indicated backflow for several of the density-unstable conditions. Generally, a larger imposed flow velocity, corresponding to a higher densimetric Froude number, is seen to suppress the tendency for backflow of the light fluid. We categorize our density-unstable displacement experiments in accordance with earlier work on buoyant pipe displacements, and find good qualitative agreement when using similar classification rules. Future work will attempt to refine the experimental methodology, and to extend the parameter regime to inclined and eccentric annuli, as well as fluid pairs with a viscosity contrast.

Although the current study is limited to the vertical and concentric annulus, this is considered highly relevant for the field application where reverse circulation cementing is known to be performed in a second stage cementing, e.g., above a weak formation [4]. Preventing backflow, i.e., ensuring instantaneous displacements should be considered an important aspect of reverse circulation job design. As such, maintaining a high imposed flow velocity and minimizing the density difference between successive fluids can aid in this regard. The previously considered field case, discussed in Refs. [4,6], suggested that it is indeed possible to achieve such conditions with conventional muds and cement slurries, and to obtain cement coverage even in the case of highly eccentric annuli. Since reverse circulation cementing remains a candidate placement strategy for geothermal wells and as a second stage cementing method for onshore or platform wells, future work should attempt to derive more precise job design guidelines to support operations and improve annular cement quality in such wells.

This work was supported by the Research Council of Norway through project number 294815.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.