A Selection Flowchart for Micromodel Experiments Based on Computational Fluid Dynamic Simulations of Surfactant Flooding in Enhanced Oil Recovery

A selection flowchart that assists, through Computational Fluid Dynamics (CFD) simulations, the design of microfluidic experiments used to distinguish the performance in Chemical Enhanced Oil Recovery (CEOR) of two surfactants with very similar values of interfacial tension (IFT) was proposed and its use demonstrated. The selection flowchart first proposes an experimental design for certain modified variables (X→: porosity, grain shape, the presence of preferential flowing channels, and injection velocity). Experiments are then performed through CFD simulations to obtain a set of response variables (Y→: recovery factor, breakthrough time, the fractal dimension of flow pattern, pressure drop, and entrapment effect). A sensitivity analysis of Y→ regarding the differences in the interfacial tension (IFT) can indicate the CFD experiments that could have more success when distinguishing between two surfactants with similar IFTs (0.037 mN/m and 0.045 mN/m). In the range of modifiable variables evaluated in this study (porosity values of 0.5 and 0.7, circular and irregular grain shape, with and without preferential flowing channel, injection velocities of 10 ft/day and 30 ft/day), the entrapment effect is the response variable that is most affected by changes in IFT. The response of the recovery factor and the breakthrough time was also significant, while the fractal dimension of the flow and the pressure drop had the lowest sensitivity to different IFTs. The experimental conditions that rendered the highest sensitivity to changes in IFT were a low porosity (0.5) and a high injection flow (30 ft/day). The response to the presence of preferential channels and the pore shape was negligible. The approach developed in this research facilitates, through CFD simulations, the study of CEOR processes with microfluidic devices. It reduces the number of experiments and increases the probability of their success.


Introduction
Oil extraction is becoming more dependent than ever on enhanced oil recovery (EOR) to improve oil production [1]. One of the best-known EOR methods is the injection of surfactants to reduce the IFT between the oil and the displacement fluid, increase the capillary number, separate the crude oil from the reservoir walls, and increase the recovery factor [2]. Surfactant injection has been a widely used technique to increase the recovery factor in hydrocarbon reservoirs. This technique has been extensively studied to optimize its application. One problem present in the surfactant injection is the interaction that it may have with a porous medium, causing adsorption of the fluid on the rock and decreasing the effective concentration of surfactant that will act at the interface between the displacing fluid and the crude. Another point of interest in the surfactant application is the change in wettability that it can cause in the medium, favoring the recovery process.
Regarding the study of surfactants, research focused on optimizing its effect using other substances that decrease adsorption in a porous medium, increasing its effectiveness in reducing IFT, or increasing the displacing fluid's viscosity. The use of ionic liquids as surfactant adsorption inhibitors in porous media have been investigated [3]. In addition, novel technologies such as nanotechnology have also been involved by applying nanomaterials of different nature (hydrophilic and hydrophobic) to avoid the adsorption of surfactants in a porous medium [4][5][6][7][8][9][10]. On the other hand, surfactants of natural origin that are more effective than synthetic surfactants and are low cost, have high availability, and are biodegradable have been investigated [11][12][13][14]. These investigations have demonstrated the great interest of the academy, industry, and the community in general focused on improving EOR processes for evaluating surfactants in search of more efficient processes.
Evaluation of the performance that surfactants may present at reservoir conditions is commonly conducted in core-flooding tests that typically do not allow for flow visualization and are of an extended duration, of the order of 80 h or more [15]. Microfluidic devices can complement, and in some cases even replace, core-flooding tests because they demand shorter evaluation times and allow for the visual characterization of the flow at the pore level [16]. Micromodels of different materials can be fabricated, such as glass [17][18][19][20][21], quartz [22], silicon [19,21], or polymers [23]. The porous media geometry in micromodels can take different forms: perfectly regular [24], partially regular [25][26][27], fractal [28], and irregular [29]. The geometry of a micromodel device can resemble a reservoir prototype that mimics the texture of scanning electron microscopy (SEM) images from a reservoir rock [30]. The most common fabrication methods for micromodels are optical lithography [31], etching [20,22], stereolithography [31], and soft lithography [32][33][34][35][36]. These fabrication procedures need to guarantee the material's transparency and aim to replicate the flow conditions. The selection of the best combination of microfluidic materials, fabrication methods, and experimental conditions can be overwhelming. There were no guides on selecting all these variables in the refereed literature for a specific microfluidic device application [28]. One way to assist this process is through a selection flowchart that indicates steps to select and classify elements according to defined criteria. Selection flowcharts are part of the decision-making processes in multiple industries such as management [37], waste treatment [38], transportation [39], equipment selection [40], and mining [41]. Moreover, in the oil and gas industry, these selection flowcharts have been used to evaluate methodology for selection of equipment used for the treatment of gas and oil [42], and to evaluate upstream water treatment [43].
This study proposes a selection flowchart that defines key geometric parameters and experimental conditions of a microfluidic test that could make the detection of differences in the performance of surfactants in the low interfacial tension range faster and more reliable. To this aim, it takes advantage of Computational Fluid Dynamics (CFD) [24][25][26][27]29,44,45] to represent the flow in microfluidic devices. In this way, this research shows the application of a selection flowchart together with CFD simulations to evaluate a set of variables considered in the experiments carried out in microfluidics to evaluate the performance of surfactants. Additionally, it provides the ability to choose a set of conditions that allows to better compare the effect of surfactants on variables of interest to the oil industry such as the oil recovery factor and other types of variables that are easy to evaluate through micromodels such as the distribution of fluids.
Although this research focuses on applying a selection flowchart, it can be used for a specific oil reservoir. In this case, the variables must be chosen within a range where Processes 2021, 9,1887 3 of 24 they are found, and a multiphase model in CFD that contemplates all the phenomenology around these variables must be used.
To engineer the selection flowchart, typical microfluidic geometries were meshed in CFD. Their performance was analyzed based on characteristic metrics for EOR processes with surfactant injection such as recovery factor, pressure drop, and breakthrough time, as well as with some that are not that frequently used but give more information on the flow and can be easily calculated in CFD, such as the fractal dimension of the flow pattern and the amount of trapped oil. Two surfactants in the low IFT range were evaluated. Sensitivity analysis was carried out to propose a flowchart that can be followed to determine the micromodel experiment that allows for better differentiation of the performance of various surfactants in the low IFT range.
A flowchart for the design of microfluidic experiments for the evaluation of surfactants, such as the one described in this study, should state the values of the inputs ( → X) which would have the most significant effect on the outputs ( → Y ), so that the performance of surfactants in the low interfacial tension range can be evaluated quickly. Figure 1 summarizes this approach. The start of the flowsheet should be a characterization of the surfactants to be analyzed based on their interfacial tension (IFT). In a second step, the user should define the modified variables of the micromodel ( → X) that will be selected in the analysis. Knowledge of ( → X) allows for an experimental design step where these variables are combined, for instance, by factorial design. In a third step, the CFD simulations of the microfluidic "experiments" proposed in the previous step are conducted. This process involves multiple CFD simulations that yield the output vector ( → Y ). A sensitivity analysis of the outputs to the IFT then provides the required information to recommend a microfluidic experiment that can best distinguish the performance of the surfactants. An extensive flowchart is shown in the Figure S10 in Supplementary Material information that includes the steps to follow when CFD simulations are not available and when evaluating one or multiple outputs ( → Y ).

Surfactants
The first step in Figure 1 involves knowledge of the IFT of the selected surfactants. The IFTs for the surfactants in the low interfacial range in these studies were: 0.037 mN/m and 0.045 mN/m. While 0.037 mN/m represents the actual IFT of a surfactant [62], a value 20% higher (0.045 mN/m) was arbitrarily selected to test the ability of the flowchart to recommend experimental conditions to differentiate the performance of both surfactants.

Modifiable Variables ⃗
If unlimited resources and time were available, all possible input and output variables could be used in the experimental design and sensitivity analysis, respectively. In fact, with the constant improvement in computational capacity and the advances in machine learning and artificial intelligence, it could be a realistic possibility in the following years. Nevertheless, in a more realistic framework, the number of variables in ⃗ and ⃗ should be defined based on the number of available resources. To assist in the process of variable selection, Table 1 proposes a scale that grades the relation between ⃗ and ⃗ . In this way, the total number of possible combinations of ⃗ and ⃗ can be reduced. In Table 1 the intensity of the greyscale is indicative of the relationship between variables. An intense gray indicates a strong relationship, while white implies almost no relation between both. For instance, the presence of a preferential flowing channel, the injection velocity, and the pore-throat connectivity have the highest effect on the oil recovery factor. In contrast, the effect of grain and pore shape is lower. The relationships shown between ⃗ and ⃗ are obtained from a literature search on the interaction of different variables in EOR processes with surfactant injection [17,[25][26][27]29,30,45,[57][58][59][60]73]. Table 1 represents a direct relationship between the modifiable variables and the outputs; for a better interpretation, it should be clarified that the combined effects of different modifiable variables must be considered. For this purpose, it is recommended to carry out tests where combinations can be made.

Surfactants
The first step in Figure 1 involves knowledge of the IFT of the selected surfactants. The IFTs for the surfactants in the low interfacial range in these studies were: 0.037 mN/m and 0.045 mN/m. While 0.037 mN/m represents the actual IFT of a surfactant [62], a value 20% higher (0.045 mN/m) was arbitrarily selected to test the ability of the flowchart to recommend experimental conditions to differentiate the performance of both surfactants.

Modifiable Variables
If unlimited resources and time were available, all possible input and output variables could be used in the experimental design and sensitivity analysis, respectively. In fact, with the constant improvement in computational capacity and the advances in machine learning and artificial intelligence, it could be a realistic possibility in the following years.
Nevertheless, in a more realistic framework, the number of variables in ( → X) and ( → Y ) should be defined based on the number of available resources. To assist in the process of variable selection, Table 1 proposes a scale that grades the relation between ( → X) and ( → Y ). In this way, the total number of possible combinations of ( → X) and ( → Y ) can be reduced. In Table 1 the intensity of the greyscale is indicative of the relationship between variables. An intense gray indicates a strong relationship, while white implies almost no relation between both. For instance, the presence of a preferential flowing channel, the injection velocity, and the pore-throat connectivity have the highest effect on the oil recovery factor. In contrast, the effect of grain and pore shape is lower. The relationships shown between ( → X) and ( → Y ) are obtained from a literature search on the interaction of different variables in EOR processes with surfactant injection [17,[25][26][27]29,30,45,[57][58][59][60]73]. Table 1 represents a direct relationship between the modifiable variables and the outputs; for a better interpretation, it should be clarified that the combined effects of different modifiable variables must be considered. For this purpose, it is recommended to carry out tests where combinations can be made.
According to Table 1, when a user wants to address the effect of the surfactant on the recovery factor, the input variables that should be analyzed are the presence of preferential flowing channels, injection velocity, and pore-throat connectivity. These are the modifiable variables that significantly impact the recovery factor (strong influence). A similar approach can be used for other outputs. Implicit in the above analysis is the a priori knowledge of ( → Y ), the vector that involves the results that are going to be measured after the microfluidic experiment. ( → Y ) depends, obviously, on the experimental setup and on the objectives of the research. Measurable variables such as recovery factor, breakthrough time, pressure drop, and fluid distribution are standard in most microfluidic experiments. In contrast, other results such as fractal dimension and displacement micromechanisms are more sophisticated and demand more time and resources from the experimentalists. Variables such as emulsion formation and drop shape respond to a particular interest in emulsification processes. While the variables Table 1 are those of more common use in the microfluidic experiments of the authors, other variables, such as species diffusion, fluid mixing, asphaltene deposition, surfactant adsorption, microemulsion properties, salinity changes, and surfactant solubility, can be readily included in Table 1. However, these have their limitations in the application of CFD models.
From the nine variables in Table 1, five were selected to illustrate the ability of the proposed flowchart to select the microfluidic experiment that distinguishes the effect of the two surfactants. (1) Recovery factor: the amount of crude oil that can be obtained from the displacement of a fluid in a porous medium, usually reported as a percentage of the initial crude oil. The recovery factor is one of the most used variables in microfluidic experiments, given that an increase in the recovery factor is the ultimate goal of any CEOR process.

Experimental Design
An experimental design should be used to define a proper combination of the variables in ( → X) to assess the applicability of the surfactant. A two-level factorial experimental design with five factors (2 5 factorial designs) was applied [74]. The two limits of each modifiable variable in the CFD simulations were selected from typical values available in the literature. A circular (−) and an irregular (+) grain shape were considered as the former makes the detachment of crude oil easier compared to a quadratic or triangular shape [29], and an irregular pore shape is a better representation of the porous media of the reservoir [30,36]. The minimum and maximum porosities were 0.5 (−) and 0.7 (+), respectively. Although these values are high compared to the typical porosity in the reservoir, much of the research in microfluidics has been conducted within this range [25,28,30,58,59,[75][76][77][78][79][80][81][82]. The minimum and maximum injection velocities were 10 ft/day (−) and 30 ft/day (+), respectively. The low-level injection velocity was chosen due to its common use in microfluidics processes [25,62,75,83]. The high level is of interest because it is indicative of the effect of high injection velocities on the flow pattern [75]. The presence (+) or absence (−) of preferential flowing channels assesses the response of other variables, such as pressure drop, recovery factor, and fractal dimension [57,59,60], to disruptions in the grain pattern. Table 2 details the 2 5 factorial experimental designs.

Sensitivity Analysis
The response in an output variable (Y i ) to changes in IFT was used as indicative of the capacity of the micromodel to differentiate the properties of the surfactant. This response was measured as a normalized sensitivity coefficient (χ i,IFT ) such as that described in Equation (1): is the maximum value of the output variable in the CFD experiment i, Y i1 is the value of the output variable in the CFD experiment i when using surfactant one, Y i2 is the response variable in the CFD experiment i when using surfactant two, and IFT 1 and IFT 2 are the interfacial tensions between each surfactant and the crude oil.

Numerical Implementation
The center in the selection flowchart in Figure 1 is the experiment based on CFD simulations (CFD experiments). In the case of this research, those CFD experiments are simulations that typically include a geometry or representation of the physical space where simulations take place, i.e., a mesh that discretizes the geometry so that the balance equations can be solved, and the actual CFD solution.

Geometry
The micromodels were considered in a two-dimensional space given their negligible depth (0.099 mm) when compared to their other dimensions (12.7a × 26.5 mm 2 ). The experimental design in Table 2 results in eight different micromodels, given the possible combinations between grain shape, porosity, and the presence of preferential flowing channels.
The geometries with circular pore shapes were generated using Matlab to randomly distribute non-overlapping circles with a radius between 0.4 mm and 0.6 mm for porosities of 0.7, and a radius between 0.5 mm and 0.7 mm for porosities of 0.5. Micromodels with irregular pore shapes were based on a micromodel template obtained from the literature [62]. A micromodel without the presence of a preferential flowing channel with an irregular pore shape is the same as that used to validate the CFD results below [62]. Preferential flowing channels were placed in the center of the porous medium with an average width of 1.8 mm positioned at an angle of 45 • with respect to the direction of the fluid. Figure 2 and Table 3 show the geometries and the main characteristics of all the micromodels used in the simulations.

Mesh
The geometry was discretized in two-dimensional, unstructured meshes with triangular elements. The evaluation of grid independence was carried out for water injection with a velocity inlet of 10 ft/day. The relative difference (δ mesh ) between the pressure drops, expressed as described by Equation (2), within micromodels with a different number of elements was considered when comparing the grids where ∆P Finer and ∆P Coarser are the pressure drops between the inlet and outlet of the microfluidic device in the finer and coarser meshes, respectively. Table 4 shows the number of cells and nodes, the pressure drop, and δ mesh for all grids. The ratio between the number of elements in the mesh i (n i ) and those in the mesh with the lowest number of elements (n 0 ) varied between 1.6 and 2.4 for an initial refinement and between 3.5 and 6.4 for final refinement. Grid independence was carried out for all eight micromodel geometries. When δ mesh was lower than 0.060, the grid with the lower number of cells was selected as the grid for simulating the surfactant injection process. Figure 3 shows the grid of micromodel (d) as an example. Figures S1-S8 in the Supplementary Material information depict the mesh for all the geometries. The geometry of the cell was adjusted to guarantee that the skewness had average values below 0.28 and maxima below 0.98 for all the micromodels.

CFD Implementation
The commercial software Ansys Fluent 19.1 was used to model the flow in the micromodels [84].
The characteristics of the machine used for the simulations are as follows:

CFD Implementation
The commercial software Ansys Fluent 19.1 was used to model the flow in the micromodels [84].
The characteristics of the machine used for the simulations are as follows: •  (3)), momentum conservation (Equation (4)), and Volume of Fluid (VOF) (Equations (5) and (6)) to represent the flow in the porous medium. The multiphase VoF approach tracks the interface between the oil and the displacing fluid. VoF calculates the volumetric fraction of each of the phases within each cell in the domain [84], taking into account the IFT between phases and the contact angle of each phase to describe the wettability of the medium: where → u = (u, v) is the velocity vector and ρ is the fluid volume-averaged density, which is considered constant due to the low compressibility of the fluids.
where p is the pressure, µ the dynamic viscosity coefficient, and → F is the vector representing external forces, which for this research is the surface tension force.
where S αi is the source term of phase i, ρ i is the density of phase i, → u i is the velocity vector for phase i, and α i is the volume fraction of phase i. The phases within the micromodel are considered immiscible. The volume fraction, α i , is 0 if the cell is empty, varies between 0 and 1 if the interface is located in the cell, and is 1 if the cell is filled with phase i. The sum of the volumetric fractions of each phase in the domain must be one.
The density (ρ) and viscosity (µ) of the fluid in a cell where the interface is located are calculated as the volume-weighted average of each of the phases (Equations (7) and (8)).
To calculate the surface tension force ( → F), the continuum surface force model (Equation (9)) is used [85]: where σ ij is the coefficient of surface tension, which represents the effect of each surfactant through the IFT, and it is considered constant and κ i is the curvature of the interface taken from inside phase i. κ i is defined in terms of the divergence of the unit vector (n) of the gradient of the volume fraction of phase i as Equations (10) and (11) describe: where n is the vector normal to the interface surface, and the phase wettability and the contact angle are used to adjust the interface curvature in areas close to the grains of the micromodel.n =n w cos θ i +t w sin θ i wheren w andt w are unit vectors that are normal and tangential to the grains, respectively, and θ i is the contact angle of phase i with the grains. From the equations presented in this section, the effect of IFT (Equation (9)), wettability of the medium (Equation (12)), phase distribution (Equation (5)), viscosity, and velocity are related, allowing the analysis of the process through solving equations in different scenarios.

Solver and Boundary Conditions
The CFD simulation used the Pressure Implicit with Splitting of Operators (PISO) approach, a convergence criterion of 0.001 for all the parameters (continuity and velocities), under-relaxation factors of 0.7, 1.0, 1.0, and 0.3 for pressure, density, body forces, and momentum, respectively, and a second-order upwind scheme.
The boundary conditions applied to the simulations were a uniform velocity inlet (left side of the geometry in Figure 2) that was varied according to the CFD experiments defined in Table 2 and a constant pressure outlet (right side in Figure 2) that was set equal to the atmospheric pressure. All grain surfaces were considered walls with total wettability to oil (θ oil = 0). The gravitational forces were supposed to have a negligible effect on the flow through the porous media. The micromodels were initially saturated with crude oil. A variable time-step approach was used with a global Courant number of two. The Supplementary Material information (Section S2) presents additional parameters of the variable time-step algorithm.

CFD Model Limitations
The chemistry of the EOR processes is fundamental because it allows the description of the processes and phenomena to understand the interaction between the phases better. Processes such as adsorption of surfactant in the porous medium, microemulsion properties, surfactant solubility, distribution of surfactant between phases, among other processes, were not considered in the development of this research. However, these are of great importance in EOR processes with surfactants, as can be seen in additional research [3][4][5][6][7][11][12][13][14]86].
This model was specified with a constant IFT (Equation (9)), and the adsorption process would modify the IFT, making it sensitive to changes in adsorption rate and fluid distribution inside the micromodel. In the model developed in this research, only the IFT effect was considered; for this reason, a constant IFT and a simple relationship with wettability were proposed.
The model developed focuses on the displacement of fluids and their interaction through interfacial tension. Parameters such as retaining surfactants in the walls can be modeled through boundary conditions and dynamic changes in IFT. For this purpose, it is essential to consider the experimental adsorption kinetics and the adsorption isotherms of the different components of the system in the medium and their effect on the effective concentration of surfactant in the displacement front [87]. Moreover, asphaltene deposition can be modeled through reaction models, where equilibrium constants are necessary [88]. Furthermore, emulsions can be modeled with an Eulerian type multiphase model for each phase present in the system [89]. However, these were not considered for the simplicity of the model and the practicality of micromodel evaluation, as done in other research [26,27,29,45].
Supplementary Material information (Section S4) gives more information about the phase behavior and the modeling carried out in this research.

Evaluation of Output Variables
The oil recovery factor was obtained through an area-weighted average of the mass of oil inside the micromodel at the end of the CFD experiment. This number was compared with the amount of crude oil at the beginning of the displacement process to compute the recovery factor. The breakthrough time was determined as when the concentration of displacing fluid at the exit of the micromodel changed from 0 to 0.03 and was reported as injected pore volume. The pressure drop was determined as the difference from line averages of the pressure at the exit and the entrance of the simulation domain. The calculation of the fractal dimension of the flow pattern at the breakthrough time involved image analysis and the fractal box-counting method [90]. With the box-counting method, the fractal dimension is obtained as the slope of the line of the logarithm of the number of boxes occupied by the pattern (N) and the logarithm of the inverse of the box size (r). Supplementary Material information (Section S3) gives more details on the procedure to compute the fractal dimension of the flow pattern. High values of the fractal dimension (close to two) indicate a more uniform displacement front, where the flow pattern looks like a square. In contrast, low values (close to one) suggest a line-like displacement front, where the fluid is not distributed throughout the available space in the micromodel. For the entrapment effect, the proportion of area enclosed by the displacing fluid was calculated. As an example, Figure 4 illustrates the entrapment effect as the oil, in red color along with the grains where the crude is stuck, is enclosed by the displacing fluid (black) on the porous media (white). The flow direction is from left to right. The fraction of oil entrapped was calculated by a custom-made subroutine that, through color differentiation, identifies areas of the micromodel where the oil is surrounded by displacing fluid and/or around a pore through which the displacing fluid has passed, the subroutine calculates the area of this section and compares it with the total area available in the micromodel, subtracting the area of the grains only to consider the amount of crude stuck to them. Thus, the areas in red color presented in Figure 4 are subtracted from the grains located in those red zones to calculate the total amount of crude oil stuck in the micromodel. The entrapment effect in Figure 4 is 0.041. in red color presented in Figure 4 are subtracted from the grains located in those red zones to calculate the total amount of crude oil stuck in the micromodel. The entrapment effect in Figure 4 is 0.041.

Fluid Properties
The properties of the displacing fluid (surfactant solution) and the oil used in the analysis were taken from reference [62]. Table 5 presents the values of density and rheological parameters for both fluids. Equation (13) shows the power-law used to estimate the non-Newtonian viscosity: = kγ (13)

Fluid Properties
The properties of the displacing fluid (surfactant solution) and the oil used in the analysis were taken from reference [62]. Table 5 presents the values of density and rheological parameters for both fluids. Equation (13) shows the power-law used to estimate the non-Newtonian viscosity: (13) where η is the apparent viscosity of the fluid, k is the consistency index that is a measure of the average viscosity of the fluid, . γ is the shear rate, and n is a measure of the deviation of the fluid from Newtonian. The properties of the oil phase were the same as those reported in [62]. The modeling of properties such as viscosity allows for obtaining a better detail of the flow characteristics within the micromodel, considering that the analyzed medium presents geometric characteristics that dispose the fluid at different stresses.

Validation of Numerical Results
The oil recovery factor, the breakthrough time, the fractal dimension of the flow pattern, and the entrapment effect at the breakthrough time were compared with experimental data reported in reference [62], and these variables were also analyzed in different scenarios in order to evaluate the expected trend with the results of the numerical simulation. As the mentioned study [62] does not report the pressure drop, the CFD predictions were compared to Darcy's law.
The experimental data were carried out in previous research [62] which described the detailed experimental setup that consisted of the following elements: a digital camera, microfluidic device, light source, computer, OEM (Original Equipment Manufacturers) syringe pump, and waste storage.
The fluids used in the experimental test consisted of synthetic brines formulated based on saltwater of a Colombian field. The brine composition consisted of 6.46 g L −1 of NaCl, 0.136 g L −1 of CaCl 2 ·2H 2 O, and 0.20 g L −1 of MgCl 2 ·2H 2 O. The surfactant used in the experimental test consisted of a mixture of hydrophilic and hydrophobic surfactant in a ratio of 80:20. In preparing the surfactant mixture, the hydrophilic surfactant was first added to the synthetic brine and then the hydrophobic surfactant was added.
Displacement tests were carried out at atmospheric pressure at a temperature of 25 • C in a micromodel made of polydimethylsiloxane (PMDS). Details of the process of fabrication can be found in [91]. The tests were evaluated through image analysis taken in a high-resolution digital camera, where the high contrast between the phases is taken advantage of, and by pixel analysis, it is possible to calculate the recovery factor and the distribution of phases within the micromodel. Other details can be consulted in [62].
To validate the predictions for the recovery factor and breakthrough time, the surfactant injection process was simulated in the oil-saturated micromodel (d) with an injection velocity of 10 ft/day, the interfacial tension between the fluid and crude oil was 0.03 mN/m; these conditions were the same as those in the experimental test in [62], as well as those of the oil and surfactant already reported in Table 5. Figure 5 shows that the numerical results present the same trend as that of the experimental data and that the change in slope distribution of phases within the micromodel. Other details can be consulted in [62] To validate the predictions for the recovery factor and breakthrough time, the su tant injection process was simulated in the oil-saturated micromodel (d) with an inje velocity of 10 ft/day, the interfacial tension between the fluid and crude oil wa mN/m; these conditions were the same as those in the experimental test in [62], as w those of the oil and surfactant already reported in Table 5. Figure 5 shows that the nu ical results present the same trend as that of the experimental data and that the chan slope around 0.55 PVI was captured. The relative errors were 10% (experiments: 0.50 ulation: 0.45) for the oil recovery factors at the breakthrough time and 12.7% (experim 0.47, simulation: 0.41) for the breakthrough time. These errors are within the uncer of the experiments.  [62]. Reproduced with permission from Céspedes-Zuluaga S, Computational fluid dynamics as a tool for the design of micromodels for the evaluation of surfactant injection in enhanced oil recovery processes; published by Universidad Nacional de Colombia, 2020. Figure 6 shows the distribution of fluids at the breakthrough time in the CFD experiments ( Figure 6a) and as predicted by the CFD simulation (Figure 6b). While the pattern is not the same, both images present a certain resemblance. The images in Figure 6 were analyzed to calculate the fractal dimension of the flow pattern as explained above. The fractal dimension for the experimental test was 1.85, while this value for the numerical simulation was 1.83 as in the experimental test, the fluid maintains a more defined front of advance with less interdigitation, while in the numerical simulation, interdigitation is more evident.  Figure 6 shows the distribution of fluids at the breakthrough time in the CFD experiments ( Figure 6a) and as predicted by the CFD simulation (Figure 6b). While the pattern is not the same, both images present a certain resemblance. The images in Figure 6 were analyzed to calculate the fractal dimension of the flow pattern as explained above. The fractal dimension for the experimental test was 1.85, while this value for the numerical simulation was 1.83 as in the experimental test, the fluid maintains a more defined front of advance with less interdigitation, while in the numerical simulation, interdigitation is more evident. The simulation response to changes in IFT, an essential part of the selection flowchart given the presence of the derivative in the sensitivity analysis, was validated with experimental data. The predicted recovery factor, breakthrough time, fractal dimension of the flow pattern, and entrapment effect were compared in a surfactant injection process in the oil-saturated micromodel (d) with an injection velocity of 10 ft/day and a value of IFT The simulation response to changes in IFT, an essential part of the selection flowchart given the presence of the derivative in the sensitivity analysis, was validated with experimental data. The predicted recovery factor, breakthrough time, fractal dimension of the flow pattern, and entrapment effect were compared in a surfactant injection process in the oil-saturated micromodel (d) with an injection velocity of 10 ft/day and a value of IFT between the fluid and crude oil of 2.7 mN/m. Figure 7 compares the distribution of fluids at the breakthrough time predicted by CFD with the one reported in the experiments in [62]. Interdigitation is evident, in the upper part of the figure, for both images. The model predicts a second preferential path at the lower part of the image, although this one is more advanced in the simulation.
The simulation response to changes in IFT, an essential part of the selection flowchart given the presence of the derivative in the sensitivity analysis, was validated with experimental data. The predicted recovery factor, breakthrough time, fractal dimension of the flow pattern, and entrapment effect were compared in a surfactant injection process in the oil-saturated micromodel (d) with an injection velocity of 10 ft/day and a value of IFT between the fluid and crude oil of 2.7 mN/m. Figure 7 compares the distribution of fluids at the breakthrough time predicted by CFD with the one reported in the experiments in [62]. Interdigitation is evident, in the upper part of the figure, for both images. The model predicts a second preferential path at the lower part of the image, although this one is more advanced in the simulation.
The simulation finely predicts the breakthrough time (0.26 PVI both for simulation and experiment) as well as the recovery factor (experiment: 0.375, simulation: 0.362), the fractal dimension of the flow pattern (experiment: 1.72, simulation: 1.76), and entrapment effect (experiment: 0.0397, simulation: 0.0417). These results give confidence in the model's performance to exhibit the effect of changes in IFT on the characteristics of the porous flow.  A final validation involved comparing the pressure drop predicted by CFD and that predicted by Darcy's law [92] for the medium's permeability (5.71 D) and the length of the porous media (26.5 mm). Both values (Darcy's law: 11.76 Pa, simulation: 11.96 Pa) are in good agreement.
On the other hand, Table 2 shows the capillary number for each system evaluated. This number indicates the relationship between viscous and surface forces and was calculated as shown in Equation (14) [93]: where C N is the capillary number, µ (N/m) and U (m/s) are the viscosity and velocity of displacement phase, σ (N/m) is the IFT between the phases inside the micromodel, and θ is the oil phase contact angle. The numerator involved in Equation (14) refers to viscous forces, while the denominator to surface forces or those related to the interface between the fluids. Due to the conditions proposed for the simulations, the capillary number takes four values: 0.013, 0.039, 0.011, and 0.032. All these values indicate the superiority of the interfacial forces (C N < 1) [65], giving importance to the interfacial tension value between the crude oil and the displacing phase. In this way, in each case evaluated, it is guaranteed that the surfactant solution is of great importance.
Considering that the decrease in IFT between the phases is one of the main mechanisms used in EOR processes with surfactant injection, as shown [2,86], this research considers this effect as the main one for the development of the selection flowchart. Table 2 presents the values of ( → Y ) for the 32 simulations proposed in the factorial experimental design. The 32 CFD experiments in Table 2 can be grouped in four cases based on the characteristics of the micromodel (grain shape and the presence of a preferential flowing channel): case 1, micromodels without the presence of a preferential flowing channel with circular grain shape; case 2, micromodels without the presence of a preferential flowing channel with irregular grain shape; case 3, micromodels with a preferential flowing channel with circular grain shape; and case 4, micromodels with a preferential flowing channel with irregular grain shape. Table 6 presents the values of the normalized sensitivity coefficient χ i,IFT (Equation (1)) obtained from the CFD results in Table 2. The last column in Table 6 Table 6 allows for identifying for each displacement condition and micromodel, the output most affected by the difference in interfacial tension presented for each surfactant. In addition, it can be identified in each micromodel which is the output where the IFT change has the most impact. This type of analysis, the result of the application of a selection flowchart, together with CFD simulations and a statistical treatment, allows for having a clearer idea regarding the necessary parameters in a micromodel for the evaluation of surfactants. A similar analysis for all the other cases can be conducted to yield Table 7, which summarizes what experimental conditions have the most significant change based on differences in IFT. This table considers the conditions where the performance of the evaluated surfactants can be better differentiated, grouped according to the grain shape and the presence of preferential flow paths. In this way, it is possible to identify parameters that can be used in microfluidic experiments for the evaluation of surfactants, highlighting some of the outputs (recovery factor, breakthrough time, fractal dimension, pressure drop, entrapment effect). For case 4 (irregular grain shape with the presence of preferential flowing channel), low injection velocity and low porosity are preferred to distinguish the effect of IFT on the recovery factor. At the same time, changes in the pressure drop become more evident for a high injection velocity and low porosity. The data in Table 7 can readily be used to determine experimental conditions that can favor the evaluation of the effect of IFT on any of the output variables ( → Y ) considered in this study. Table 8 presents a global summary of the results in Table 7 that provide a comprehensive response to the question of what experimental configuration would best differentiate between two IFTs. While Table 7 responds to the question of what conditions are preferred to have a more noticeable effect on an individual output variable (Y i ), Table 8 considers the set of experimental conditions that would indicate the largest change in ( → Y ) for smallest differences in IFT. Table 8 indicates that for case 3, for instance, a high injection velocity and low porosity, in the range of the simulations used in this paper, should be preferred when evaluating two different surfactants.

Application of the Micromodel Experiment Selection Flowchart
Finally, from the values of χ TOTAL ,IFT in Table 6, it is possible to choose the micromodel among the eight evaluated in this research. The conditions where a higher value of χ TOTAL ,IFT (76.98) was obtained were those of the micromodel with circular grain shape, with the presence of a preferential flowing channel, a porosity of 0.5, and an injection velocity of 30 ft/day. This means that this micromodel with the evaluated conditions is the one that allows the best differentiation of the performance between the surfactants analyzed in this research. This was obtained from the application of the selection flowchart. The example outlined above illustrates how the flowchart proposed in Figure 1, which uses CFD simulations to represent the flow in a microfluidic device, can assist in selecting the microfluidic experiment that would have more success in distinguishing the performance of two surfactants with very similar IFTs. Also, the conditions where the micromodel best differentiates the performance of surfactants is found for the highest capillary numbers, with values of 0.039 and 0.032, corresponding to an IFT of 0.037 mN/m and 0.045 mN/m, respectively. However, this indicates a significant influence of viscous forces, specifically the injection velocity, other parameters must also be considered as will be discussed later depending on the output variables evaluated.
Further information can be obtained from the analysis of Table 6, particularly on the importance of the output variables on the evaluation of surfactants. The sum of all χ i,IFT over all the CFD experiments indicates that the entrapment effect (∑ χ Ent,IFT = 218.9) is the variable that can more easily differentiate the performance of both surfactants. The recovery factor and breakthrough time (∑ χ RF,IFT = 50.9, ∑ χ BT T,IFT = 56.6) can also be used to distinguish the performance of surfactants with very similar IFTs. The fractal dimension (∑ χ FD,IFT = 12.6) and pressure drop (∑ χ P,IFT = 11.5) have the smallest sensitivities to changes in IFT. The strong influence of the surfactant on the entrapment effect, the recovery factor, and the breakthrough time is probably because surfactants have a substantial effect on entrapment as they modify the surface interaction between the oil and the porous media, and increase the amount of crude oil displaced through the micromodel to increase the recovery factor and reduce the breakthrough time. The low effect on the fractal dimension of flow indicates that in the small length scale of a microfluidic device, such as those used in this study, the fractal dimension of flow differences are slight for surfactants of very similar IFTs, such as those used in this study. The same is true for the pressure drop.
A similar attempt at correlating the microfluidic grain shape (circular vs. irregular) or the presence or absence of preferential flowing channel does not allow a definitive conclusion as similar values of ∑ χ i,IFT were obtained for these variables    no pre f erential channel  increase the ability of the microfluidic experiment to distinguish between two surfactants. This indicates that more extreme conditions, lower porosity, and a higher velocity make the effect of the surfactant in the oil recovery process more evident.

Conclusions
This research proposed and demonstrated a flowchart for selecting microfluidic experiments to evaluate the performance of surfactants in the low interfacial tension range. The flowchart guides the user through comprehensive CFD simulations, experimental design, and sensitivity analysis with the ultimate goal of defining the microfluidic characteristics and the experimental conditions that can make the difference between two surfactants either in one specific output variable, e.g., recovery factor, or in a set of performance variables more evident.
For the output variables that can be measured in microfluidic experiments and that were considered in the analysis, i.e., recovery factor, breakthrough time, fractal dimension, pressure drop, and entrapment effect, the latter is the one that makes the differences between the surfactants used in the study more evident; however, the recovery factor and the breakthrough time can also be used to distinguish between surfactants.
For the microfluidic design characteristics and experimental conditions evaluated in this study, a high injection velocity and a low porosity have the most significant effect on the ability of the CFD experiments to differentiate the performance of the surfactant. The presence of preferential flowing channels on whether circular or irregular grains form the microfluidic does not have a significant effect when the objective is to analyze the performance of a surfactant.
Although this type of process facilitates the evaluation of EOR technologies with surfactant injection, it is recommended to carry out experimental tests to complement the results obtained from this type of research. The objective of this research is not the design of EOR processes, but of micromodels that help to select surfactants and other variables before the application of the technologies in the field.
Phenomena such as adsorption of surfactant in the porous medium, microemulsion properties, surfactant solubility, distribution of surfactant between phases, among other processes, were not considered in the model proposed. However, they are essential for evaluating surfactants and should be included in future studies.
The results obtained in this research facilitate the design of micromodels to evaluate surfactants in EOR processes. The application of surfactants in the field involves different types of phenomena and interactions between surfactant-porous media that were not considered in the modeling present in this investigation. For this reason, the results and conclusions obtained in this investigation are limited to the design of micromodels in the evaluation of surfactants.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/pr9111887/s1, Figure S1: Detail of the mesh for micromodel (a), Figure S2: Detail of the mesh for micromodel (b), Figure S3: Detail of the mesh for micromodel (c), Figure S4: Detail of the mesh for micromodel (d), Figure S5: Detail of the mesh for micromodel (e), Figure S6: Detail of the mesh for micromodel (f), Figure S7: Detail of the mesh for micromodel (g), Figure S8: Detail of the mesh for micromodel (h), Figure S10: Extensive Selection Flowchart for micromodel experiments. Table S1: Parameters for the variable time-step algorithm used in the simulations, Section S3: Fractal dimension of the flow pattern, Section S4: Phase behavior.