Automated determination of interfacial tension and contact angle using computer vision for oil field applications

Contact angle and surface tension are the two most widely used surface analysis approaches for reservoir fluid characterization in petroleum industries. The pendant drop method has among the most widely used techniques for the estimation of surface tension. The present work utilizes a python-based computer program to automatically determine interfacial tension (IFT) and contact angle from the pendant drop image acquired from a typical pendant drop apparatus. The proposed program uses python-based image processing libraries for the analysis of the pendant drop image. Also, the program is tested on images acquired from the standard solutions for the IFT and contact angle calculation showing promising results with a standard deviation of less than 1.7 mN/m.


Introduction
Surface analysis is an essential aspect of the petroleum and process industries. This helps in designing an efficient and economical process. The recovery of oil from the reservoir depends on the wettability, and the surface tension properties of the rock and crude oil contact angle and surface tension are two surface analysis approaches that can be used to characterize the reservoir fluid (Saxena et al. 2019). Hydrocarbon reservoir comprises immiscible phases within its porous structure, predominantly oil, gas, and water. Typically, there is a formation of a thin film between the two fluid surfaces. It is dependent upon the chemical structure, flow profiles, and properties of the different phases. With the dissimilar forces applied by two immiscible fluid phases at the interface, the interfacial tension (IFT) develops. This causes resistance in the miscibility of two fluids, making it difficult for the fluid to flow in a porous medium. IFT and contact angle are important parameters dictating oil recovery from petroleum reservoirs (Sheng 2013;Alnoush et al. 2019). Contact angle helps in studying the wettability of the reservoir rock. Wettability is defined as the tendency of the fluid to adhere or spread over the surface of the rock (Kesarwani et al. 2021a). It helps in characterizing the rocks as water wet or oil wet, whereas IFT is the force acting between the two interfaces of the immiscible fluid. The IFT dictates the ease of oil movement from the porous structure. The lower the IFT, the easier the hydrocarbon can be displaced from the reservoir rock's pore spaces, overcoming the capillary forces. The capillary forces keep the oil trapped owing to its more considerable interfacial tension and lower capillary diameter (Arashiro and Demarquette 1999).
Rheological theories have been developed to infer IFT from small-amplitude oscillatory shear measurements. They involve evaluating a profile of either a sessile drop, a spinning drop, or a pendant drop. Dynamic methods based on the breaking thread and embedded fiber were also used to determine IFT (Zamora et al. 2018).
The pendant drop method is one of the most advanced and widely accepted techniques for estimating IFT. The estimation of the surface tension by analyzing the shape of the drop has been the prime method. However, the results of this method are subjected to the assessment of images of the drop. The recent advancements in the field of image processing have made the job convenient. With the python-based image processing techniques, the images are converted to an appropriate color format, and the image contours' equations can be solved.

3
The current work presents a python-based program for an automated determination of the interfacial tension by processing the image of pendant drop (captured by the camera) using the library OpenCV (computer vision) (Rakesh et al. 2004). Due to this interfacial tension between immiscible fluids (one fluid tends to adhere to the pore surface over the other fluid), the reservoir rocks show a characteristic called wettability. The wettability of a reservoir rock affects the flow of reservoir fluid over its surface. The wetting phase tends to occupy the more diminutive and remain under dynamic reservoir conditions. However, water-wet rocks are desirable for higher oil recovery (Khaksar Manshad et al. 2016).
The contact angle measurement is a method to quantify the spreading capability of the fluid over the rock surface. It is defined as an angle formed by the interaction of a fluid-solid interface by applying a tangent line from the contact point along with the oil-water interface, as shown in Fig. 1. It is observed that the contact angle differs significantly depending on the wettability of the rock. When the rock's surface is water wet, water droplets tend to spread out to the boundaries due to stronger attraction between the water phase and the solid surface, leading to a contact angle significantly smaller than 90°. On the other hand, an oil-wet rock tends to have a contact angle larger than 90°. Rocks having intermediate wettability measurements are found to have a contact angle around 90° (Kesarwani et al. 2021b). The current work proposes a method for determining the contact angle using image processing techniques (Wang and Zhang 2008). The novel algorithm can estimate the contact angle using the lower-resolution images.

Literature review
Andreas et al. proposed a straightforward method for obtaining this amount in the 1940s by dividing the greatest drop diameter (D e ) by the drop diameter (D s ) determined at a distance (D e ) from the apex (Andreas et al. 1938). The Bond number may then be determined by comparing the ratio S = D s /D e to tables, yielding the interfacial tension. Andreas et al. obtained these tables experimentally. Later, these values were revised by performing numerical integration of the Young-Laplace equation (Fordham and A 1948). This method provided a straightforward way for computing interfacial tension, but it disregarded a considerable percentage of data. In 1983, two seminal publications were published that built computer methods to use all available data, considerably improving the method's accuracy (Huh and Reed 1983;Rotenberg et al. 1983). These methods compared the overall drop profile to the theoretical drop profile by calculating the sum of the squared residuals between each experimental data point and the theoretical drop profile. While the approaches are similar in many ways, Huh and Reed used an estimated expression rather than the precise equation used previously (Huh and Reed 1983;Rotenberg et al. 1983). Earlier technique treated the apex position as an indefinite quantity that was computed concurrently with radius of the apex and the shape parameter, which is same as the Bond number, resulting in higher accuracy (Rotenberg et al. 1983). Hoorfar and Neumann described the developments in computational approaches of the ADSA (axisymmetric drop shape analysis) algorithms (Hoorfar and W. Neumann 2006). Jennings and Pallas further improved on earlier technique by using rotational discrimination (Jennings and Pallas 1988). This can be viewed as a modified Gauss-Newton approach, to accomplish the optimization procedure, substantially lowering the computing time. This paper also included a rigorous error analysis that offered precise, conservative ranges for the experimental error. Despite the fact that Gauss-Newton optimization is often used and may achieve extremely fast convergence between the experimental and fitted droplet profiles, it can fail to converge if the theoretical profile's initial estimate is not reliable. The pendant drop technique has recently been improved to accommodate a larger range of drop alignments, including capillary bridges between two parallel plates (Cabezas et al. 2006;Kalantarian et al. 2011;Vega et al. 2014). These arrangements have given pendant drop tensiometry a major boost, allowing it to properly detect interfacial tension even if Bond Number reaches up to zero (Neeson et al. 2014).
This article provides description of the python program, which takes the use of certain open-source libraries to calculate the IFT from the pendant drop image. The program first determines the length and breadth of the pendant drop from the acquired image. This is done using the contour function of OpenCV library to draw a bounding box around the pendant drop image. The length of the pendant drop is then taken as D e (equatorial diameter). The image is then cropped in such a way that the part of pendant drop, which is present below the point present at the distance equal to D e from the apex of the pendant drop, is removed. The cropped image is then used to calculate D s (syringe diameter). Both of these obtained diameter (syringe diameter and equatorial Fig. 1 The schematic illustrating contact angle (θ) 1 3 diameter) were then converted to cm unit from the pixel unit. This is done using the pixel-to-cm ratio obtained using the ratio of the known diameter of the syringe of the equipment and obtained D s from the pendant drop image. Then, the ratio of D s and D e (S = D s /D e ) is used to calculate the H using the equations derived by Marvin D. Misak (Eqs. 3a-3f). This H is finally used, along with the density obtained from the user, to finally get the IFT value of the pendant drop using Eq. 1.

Theory of the pendant drop
In 1882, Bashforth and Adams proposed a new equation to derive a pendant drop's theoretical form and calculated tables of drop contours. The IFT was estimated by comparing the experimental data with the curve generated using the theoretical data. Images of the evolving drop are taken as a function of time to infer the experimental data better. However, this procedure is very complex and time taking. The process was simplified using the following empirical relationship proposed in the work by Berry (Bashforth S and Adams 2011;Berry et al. 2015): where is the interfacial tension, Δ is the density difference between the fluid phases, D e is the equatorial diameter of the drop (Fig. 2), H is a correction factor which is related to the shape factor of the pendant drop, S, defined as: where D s is the drop diameter measured horizontally at a distance D e away from the apex of the drop. ( Marvin D. Misak has derived six empirical equations (Eqs. 3a-3f) for generating the 1/H vs. S curve. The following shape factor helps in estimating the correlation parameter 1/H. The equations help in estimating the IFT using the pendant drop technique. The equations are valid for S values from 0.30 to 1.00. (Drelich et al. 2002).
for The computed 1/H values illustrate a minimal divergence (less than ± 0.00051) from the previously published data. In the usually encountered range of S values, the deviation is less than ± 0.00006 (Drelich et al. 2002).

Data acquisition and description
The pendant drop method was used for the surface tension measurement. A syringe pump capable of providing 10,000 psi supplied from DCAM Engineering, India, was connected to the 1/8-inch flowline with a syringe of known outer diameter attached at the other end. The syringe pump and a glass chamber were placed to have a light source at the back end, and the camera was attached to the front end to capture the image (Fig. 3). The idea of the experimental setup was to flow the standard solution using a syringe pump at a minimum flow rate so the liquids form a drop that could be captured and analyzed. The syringe pump was filled with the standard solutions of known concentration, and a flow rate of 0.001 ml/min was maintained in the syringe pump (Kesarwani et al. 2021c). Upon exiting the syringe, the liquid forms a pendant drop. At the time, when the contour of the drop stopped evolving any further, the image of the pendant drop was captured with the camera that is attached at the front end. The flowlines were cleaned properly by flushing distilled water through it. The complete experiment was performed at ambient conditions. The standard system used for procuring the images is given in Table 1. Figure 4 depicts the general flow of the proposed algorithm to calculate the IFT value from the acquired data.

Preparing the image for image processing
Throughout image processing, the library OpenCV (computer vision) is used to manipulate the image without changing its properties (Rakesh et al. 2004). The image of the pendant drop, which is to be analyzed, is imported into the program. It is challenging to extract the information from  The colored images are better for the visual appeal of humans, but machines focus on the applications of image processing. So, it is crucial to convert the image to grayscale as it allows the program to work with much fewer pixels than the original image while maintaining the structural aspect of the image. Grayscale is a range of shades of gray without any actual color. The darkest possible shade is black, and the lightest possible shade is white. The shades of gray in between white and black are represented by equal brightness intensities of the three primary colors (red, green, and blue). Figure 5a shows the grayscale image of the pendant drop. An edge detection function, called Canny edge detection (Morse 2000), is used on the grayscale image of the drop to obtain the binary edge of the drop. Before finding contours, converting the image into binary by applying threshold or canny edge detection is necessary (Rakesh et al. 2004). Figure 5b shows the detected edges of the pendant drop by the Canny edge detection function.

Finding contours using OpenCV and drawing bounding rectangle
Contours are very helpful in shape analysis, defining the size of the object of concern, and object detection. Contours are simply defined as the lines that join all those points along the boundary of an image with the same intensity. In OpenCV, finding contours is like finding an object which is white from a black background. Thus, the object should be found in white, and the background should be colored black. This is the reason we apply threshold or edge detection before finding the contours.
The function cv2.contourArea() is used to remove contours that do not have a large enough area. OpenCV contours have many features (Bradski et al. 2009), one of which is called the bounding box. This contour feature can be used to make a bounding box, of minimum area, around the contours fitted in the detected edges of the pendant drop image. Figure 6 shows the bounding box on the pendant drop image. The width of this box is directly used as the equatorial diameter, D e (refer Fig. 2), of the pendant drop. This D e (equatorial diameter) is used in the calculation of the IFT, σ (refer Eq. (1)), and the shape factor S (refer Eq. (2)). In the case of the pendant drop shown in Fig. 6, the width of the bounding box (219.0 px) is taken to be the equatorial diameter (D e ) for calculation of the IFT.

Cropping the image and again fitting contours to find drop diameter D s and syringe diameter
The image of the detected edge of the pendant drop is cropped vertically using OpenCV such that the cropped image has the section of the image which starts at a distance equal D e (equatorial diameter) away from the apex of the pendant drop as shown in Fig. 7.
The image is cropped such that the distance between the bottommost points of both the edges (points C and D in Fig. 7) is equal to the drop diameter D s , and the distance between topmost points of both the edges (points A and B in Fig. 7) is equal to the syringe diameter (refer "Theory of the pendant drop" section and Fig. 2).
The function that finds contour(), described in "Finding contours using OpenCV and drawing bounding rectangle" section, is used on this cropped image (shown in Fig. 7) to determine the coordinates of these four points A, B, C, and D, as shown in Fig. 8. These coordinates are then used to calculate the AB and CD lengths, which denote syringe diameter and drop diameter, respectively (refer "Theory of the pendant drop" section).

Calculating IFT
A problem that arises while making measurements from the image using programming is that the unit of length is measured in pixels (px). However, to calculate, we need a length of D e (equatorial diameter) in cm to put in Eq. (1). This   problem can be solved by using a reference length that is known in both units. As the syringe diameter is known to us or can be measured, the ratio of length AB and syringe diameter can be used as a reference scale to convert the units of equatorial diameter, D e , from px (pixels) to the centimeter.

Preparing the image for image processing
The captured image is imported and converted to grayscale. Converting to grayscale, it is important because it allows the program to work with much fewer pixels than the original image while maintaining the structural aspect of the image (refer "Preparing the image for image processing" section). An example of a grayscale image of the acquired drop image is shown in Fig. 9.
Although the grayscale image is much better than the original mage, applying threshold helps eliminate the noise and reflect the image. Thresholding creates a binary image, revealing clearly defined image objects and boundaries as either white or black. If the pixel value is greater than a threshold value, it is assigned one value (maybe white), else it is assigned another value (maybe black). Figure 10 shows an example of a thresholded image of the drop. The thresholding function uses three arguments. The first argument is the source image, which should be grayscale. The second argument is the threshold value, which is used to classify the pixel values. The third argument is the max Val, which represents the value to be given if the pixel value is more than (sometimes less than) the threshold value (Morse 2000;Rakesh et al. 2004).

Taking user inputs for measuring the contact angle
The thresholded image is plotted on a graph and shown in a window to the user in Fig. 10. The user inputs play a significant role in determining the contact angle from the image. The program asks the user to choose the two contacts points on the boundary of the drop. An example of user inputs is shown in Fig. 11 by the red plus sign ( +). The function records these two points x and y coordinates in the form (x left, y left) and (x right, y right). These points are used to determine the tangent to the drop whose angle with the horizontal gives the contact angle of the drop (Njobuenwu 2007).
Here the data are called (x left, y left) and (x right, y right). The equation gives the slope of the line connecting these two points: where m is the gradient. Given m, the contact angle of the drop can be easily calculated by: where θ is the contact angle (refer Fig. 1).

Results
Python program was used on some images of the standard liquids with the known value of IFT and contact angle. The value predicted by the program is reported in Table 1 against the actual values of the drop's IFT and contact angle. The standard deviation in each predicted IFT and contact angle value is reported to estimate the accuracy of the proposed algorithm.
Considering the standard deviations, it is clear that our program is doing a decent job of estimating the IFT value of the fluid. However, in the measurement of contact angle the measurement requires manual interference for selecting the boundary of the droplet. The standard deviation is very less showing a maximum deviation of 1.7, which is within the acceptable range for human error while performing the experiments. This indicates that the python program presented in this study is entirely consistent and efficient.
The presented algorithm for IFT estimation consists of four steps: 1) Acquisition and preparing the pendant drop image for processing. 2) Processing the acquired image to estimate different geometrical parameters required to calculate IFT. 3) Taking user input of density difference between fluid under evaluation and its surrounding. 4) Calculating and giving out the determined value of IFT.
The proposed python program provides a straightforward and reliable approach for determining surface tensions and contact angle of various substances for reservoir characterization. It has utility in data interpretation and statistical analysis, both of which are crucial in modern science. In terms of speed and accuracy, the approach utilized here is very effective and efficient.

Conclusion
The pendant drop approach for evaluating IFT for reservoir fluid characterization was discussed in our present work. The IFT is estimated from the geometric profile of a pendant drop of the liquid at mechanical equilibrium. The previous techniques involved costly equipments, still yielding in error (5) = −tan −1 (m) outside the permissible limits. These techniques were very time-consuming and inaccurate. The presented program uses the image acquired by the equipment for much more accurate and faster estimation of the IFT. The program further can be used for estimating the contact angle of the fluid over different surfaces.
Funding No external source of funding was used for this project.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.