Suppressing Viscous Fingering in Porous Media with Wetting Gradient

The viscous fingering phenomenon often occurs when a low-viscosity fluid displaces a high-viscosity fluid in a homogeneous porous media, which is an undesirable displacement process in many engineering applications. The influence of wetting gradient on this process has been studied over a wide range of capillary numbers (7.5 × 10−6 to 1.8 × 10−4), viscosity ratios (0.0025 to 0.04), and porosities (0.48 to 0.68), employing the lattice Boltzmann method. Our results demonstrate that the flow front stability can be improved by the gradual increase in wettability of the porous media. When the capillary number is less than 3.5 × 10−5, the viscous fingering can be successfully suppressed and the transition from unstable to stable displacement can be achieved by the wetting gradient. Moreover, under the conditions of high viscosity ratio (M > 0.01) and large porosity (Φ > 0.58), wetting gradient improves the stability of the flow front more significantly.


Introduction
The displacement of immiscible two-phase fluids in porous media is an important subject in natural and engineering fields, including oil exploration [1], carbon dioxide storage [2], fuel cells [3,4], batteries [5], water purification [6], general electrochemical energy storage [7,8] etc. For example, in oil exploration, the preferential flow of displacement fluids will bypass larger areas of oil, trapping some of the oil in the reservoir, thus reducing oil recovery efficiency [9]. In terms of carbon dioxide storage, supercritical carbon dioxide is injected into a deep saline aquifer. The viscous stability of the primary drainage process is of major interest for the injection of carbon dioxide in saline aquifers, since it determines the spread of the carbon dioxide plume in the target aquifer and consequently the initial utilization of the pore space for carbon dioxide storage [10]. It is of great significance for understanding the flow laws of immiscible fluids in porous media.
In order to better study the flow law of fluid in porous media, it is necessary to quantitatively analyze the results. There are several parameters for quantitative analysis, including the displacement efficiency, saturation, interface length and fractal dimension. The displacement efficiency refers to the percentage of the volume of the displaced fluid flowing out and the initial filling volume [11]. The saturation refers to the percentage of the displacement fluid volume in the total pore volume [12]. The saturation can be the percentage at different times, which is different from displacement efficiency. The length of the interface is defined as the length of the interface between the displacing fluid and the displaced fluid. It is often employed to characterize the instability of the fluid front. The longer the interface length, the more unstable the flow along the front of the fluid. The fractal dimension is described as a measure of the space-filling capacity of a pattern, which is used as an indicator to represent the complex geometric form, to compare the changes of details in the pattern, and to reflect the effectiveness of complex shapes occupying space. It is a measure of the irregularity of complex shapes. The more complex the fractal dimension is, the more complex the shape is [13].
The flow of fluids in porous media can be divided into three types: stable displacement, capillary fingering and viscous fingering [13]. There exist two dimensionless parameters by which to characterize the flow of fluids in porous media, including the viscosity ratio M = µ 1 /µ 2 and the capillary number Ca = µ 1 U/σ; here, µ 1 is the viscosity of the displacing fluid. µ 2 is the viscosity of the displaced fluid, U is the rate at which the displaced fluid is injected, and σ is the surface tension. The phase diagrams of M and Ca have been summarized [14], which is of great significance for further understanding the flow laws of fluids in porous media. Zhang [15] and Zheng [16] have also summarized the relevant phase diagrams when they study the flow of fluids in porous media. The phase diagrams summarized by these researchers are not identical in the boundary ranges of different flow patterns, which may be due to the wettability, gravity and surface roughness [17].
Zhao et al. [18] studied the influence of wetting angle between fluid and porous media, and found that the influence of wettability is not monotonic. Jung et al. [19] studied the influence of wettability by combining experiments and simulations. Their research shows that wettability has a great impact on the displacement process of immiscible fluid. When the wetting angle between the displacement fluid and the porous medium is less than 80 • , the saturation of the final displacement fluid is higher because the interface of the fluid front moves smoothly and the displacement fluid is not trapped. Even when studying the influence of wettability, the law of influence on wettability is not completely consistent. On the one hand, the pore channels in porous media are uneven, especially in the displacement problem in porous media in the engineering field, where the pores are very irregular. Even in the displacement problem in homogeneous porous media, the nonuniformity of the flow front interface will also have a great impact [20]. Golmohammadi et al. [21] used experimental methods to study the comprehensive effect of gravity and wettability on fluid flow in two-dimensional porous media, and found that gravity plays an important role in improving the filling effect. Hu et al. [22] studied the influence of roughness on the flow pattern in rock fractures and quantized the energy dissipated in the process of fluid invasion dominated by capillary force. Moreover, scholars have also made some progress in improving the stability of the fluid front flow in porous media. Rabbani et al. [11] proposed a method of suppressing viscous fingering by designing the gradual and monotonic variation of pore sizes along the front path. Lu et al. [23] investigated the influence of pore size gradient and pore-scale disorder on the displacement process when a non-wetting fluid displaces a wetting fluid, and found that a sufficiently large gradient can completely suppress capillary fingering.
With the development of science and technology, the investigation of immiscible fluid flow in porous media is not limited to experiments, and many scholars have also tried to use numerical simulation methods to carry out corresponding research. The lattice Boltzmann method (LBM) is a mesoscopic numerical simulation method between macro and micro. Compared with macro methods, it is more convenient for dealing with a complex boundary and, at the same time, it overcomes the limitations of the micro calculation method on size. It has the advantages of the automatic capture of a two-phase interface without manual processing or flexible boundary processing, and is suitable for parallel computing [24]. The lattice Boltzmann method has been widely used for simulating multiphase flow [25], micro-nano scale flow [26], turbulence [27], flow-induced vibrations [28], heat and mass transfer [29], and porous media flow problems [30,31]. Liu et al. [30] simulated the immiscible flow of wetting fluid and non-wetting fluid in two porous medias. Their simulation results confirmed that three different displacement modes are related to capillary number, viscosity ratio and the heterogeneity of porous media. Shi et al. [31] investigated the basic physical mechanism of Newtonian fluid replacing non-Newtonian fluid in porous media, which revealed the displacement mechanism of Newtonian fluid to non-Newtonian fluid from a mesoscopic perspective. Lautenschlaeger et al. [32] proposed a homogenization method for simulating multiphase flow in heterogeneous porous media, which is based on the lattice Boltzmann method and combines the gray level with the multi-component Shan-Chen method. This method makes the fluid-fluid and solid-fluid interactions in pores less than numerical discretization, and has been successfully applied to solving various single-phase and two-phase flow problems.
Unstable flows will create adverse effects for oil production, carbon dioxide storage and other aspects. In order to reveal the law of the unstable flow of fluid in porous media, some scholars have obtained the relevant phase diagrams through a large number of experiments. Some scholars have tried to improve the unstable flow of fluid by designing the structure of porous media. However, how a wetting gradient impacts viscous fingering remains unknown. In this paper, we used the lattice Boltzmann method to simulate the displacement process of fluids in porous media with a wetting gradient. The influence of the wetting gradient was studied over a wide range of capillary numbers (7.5 × 10 −6 to 1.8 × 10 −4 ), viscosity ratios (0.0025 to 0.04) and porosity (0.48 to 0.68). The flow results were quantitatively analyzed by means of quantitative parameters such as fractal dimension, displacement efficiency, saturation and interface length. Moreover, a phase diagram of flow stability results related to wetting gradient and capillary number was drawn.

Mathematical Model
The lattice Boltzmann method is used for two-dimensional numerical simulation to simulate the displacement process of fluid in porous media. The lattice Boltzmann method is a mesoscopic simulation method between micro-molecular dynamics and macrofluid dynamics. This method uses the Boltzmann transport equation to calculate two processes of collision and migration between micro-particles, replacing macro-particles with micro-particles. The macroscopic parameters of the system can be obtained by statistical averaging of a large number of particles without concern for the motion state of each particle. The motion behavior of the whole fluid is simulated and the corresponding macroscopic phenomena are analyzed. The Boltzmann transport equation is: (1) is the distribution equation of particles in the direction i of position r at time t.
Macro density ρ and velocity v on the node can be obtained by integration: The discrete velocities model adopted was the D2Q9 model, and its discrete velocities are shown in Figure 1.
The D2Q9 model illustrated in Figure 1 involves nine velocity vectors defined by:  The D2Q9 model illustrated in Figure 1 involves nine velocity vectors defined The evolution equation of the particle distribution function of the multiple-rela time model is:     Hence, the transformation matrix M is The evolution equation of the particle distribution function of the multiple-relaxation time model is: where where ρ is the density, e is the kinetic energy, is related to the kinetic energy square, j x and j y correspond to components of momentum, q x and q y correspond to the energy components, and p xx and p yy correspond to the symmetric traceless viscous stress tensors.
Hence, the transformation matrix M is S ij = diag(s 0 , s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 ) is the diagonal relaxation matrix. In this approach, the kinematic viscosity υ and the bulk viscosity µ are related to the following relaxation parameters: c s = c/ √ 3, c s is the sound velocity of the lattice; c = δx/δt, c is the lattice velocity; δx is the length step; and δt is the time step.
Applying the transformation matrix to the equilibrium probability distribution function, the raw moments at the equilibrium are The moments ρ eq , j eq x and j eq y are not required as they will be multiplied by s 0 , s 3 , and s 5 , which are zero.
The bounce-back boundary condition is used to implement the no-slip condition on a geometry wall. The incoming probability distribution functions at the wall node are reflected back to the initial fluid node, giving value to the unknown probability distribution functions. In the illustration in Figure 2, at a time t, the distribution function f 7 is reflected back to f 5 , f 3 to f 1 , and f 6 provides f 8 value.  The inlet and outlet boundary conditions used in the current simulations are ve and pressure boundary conditions, respectively. The approach consists of the formu of a linear system with mass and momentum conservation. This linear system will pr the value of the unknown probability distribution functions created after the strea step as well as the value of  if a velocity boundary condition is imposed, or u if is a pressure condition at the boundary. According to the Figure 2 example, after rear ing the moments in Equations (2) and (3) (density, x-momentum, y-momentum), the system can be written as: The inlet and outlet boundary conditions used in the current simulations are velocity and pressure boundary conditions, respectively. The approach consists of the formulation of a linear system with mass and momentum conservation. This linear system will provide the value of the unknown probability distribution functions created after the streaming step as well as the value of ρ if a velocity boundary condition is imposed, or u if there is a pressure condition at the boundary. According to the Figure 2 example, after rearranging the moments in Equations (2) and (3) (density, x-momentum, y-momentum), the linear system can be written as: A fourth equation is required to close the system and compute the value of velocity or pressure. A valid assumption is to apply the bounce-back rule for the nonequilibrium part of the particle distribution normal to the boundary. Therefore, considering that an inlet boundary is imposed on the left nodes of Figure 2, the fourth equation is Solving the linear system, the unknown probability distribution functions are For a velocity boundary condition, the linear system gives: For a pressure boundary condition, p = ρc 2 s ; taking v = 0, the u is defined by:

Simulation Setup
We used the computational fluid dynamics software XFlow to carry out the relevant numerical simulations. The phase field algorithm was employed for calculating the multiphase flow. The temperature type was set as isothermal. The computational domain and boundary conditions are shown in Figure 3. The size of the computational domain was 16 mm × 10 mm. The cylinders represent an impermeable solid material, while the area formed between the cylinders are the pore channels. The diameter of the cylinder and the spacing between cylinders are described in Figure 3A. The displaced fluid initially filled in the pores of the porous media, and the displaced fluid flowed in at a certain speed from the bottom, forcing the displaced fluid to flow out at the top. The initial gauge pressure field was 0 Pa (the corresponding actual pressure was 1.013 × 10 5 Pa) and the initial velocity field was 0 m/s. The bounce-back boundary condition was used to implement the no-slip condition on the solid walls. The inlet and outlet boundary conditions used in the current numerical simulations were velocity and pressure boundary conditions. The displaced fluid was silicon oil (the wetting phase) and the displacing fluid was water (the non-wetting phase). The parameters of these two fluids were as follows: the density and viscosity of displacing fluid were set as ρ 1 = 998.3 kg/m 3 , µ 1 = 1 mPa s. The density and viscosity of the displaced fluid were set to ρ 2 = 960 kg/m 3 , µ 2 = 200 mPa s, and the viscosity ratio was M = 0.005. The surface tension σ = 28.2 mN/m. The inlet velocity was 0.001 m/s unless otherwise specified. The wetting gradient was set as follows: Take the wetting gradient ∆θ = 4 • as an example, as shown in Figure 3B. The wetting angle of the lowest row of cylinders was set to 90 • , and the same row of cylinders had the same wetting angle. The wetting angle was increased by 4 • for each row up to the last row, thus creating the porous media of wetting gradient ∆θ = 4 • . The wetting gradients ∆θ = 1 • , ∆θ = 2 • and ∆θ = 3 • were arranged in the same way. The wetting gradient ∆θ = 0 • meant that the overall wettability of the porous media was uniform and the overall wetting angle was 150 • . otherwise specified. The wetting gradient was set as follows: Take the wetting gradient Δθ = 4° as an example, as shown in Figure 3B. The wetting angle of the lowest row of cylinders was set to 90°, and the same row of cylinders had the same wetting angle. The wetting angle was increased by 4° for each row up to the last row, thus creating the porous media of wetting gradient Δθ = 4°. The wetting gradients Δθ = 1°, Δθ = 2° and Δθ = 3° were arranged in the same way. The wetting gradient Δθ = 0° meant that the overall wettability of the porous media was uniform and the overall wetting angle was 150°. The red area is the displacing fluid, which is injected at a certain flow rate from below. The blue area is the displaced fluid, flowing out from the upper outlet. The outlet gauge pressure is 0 Pa. The horizontal line between the red area and the blue area represents the interface between the displacing fluid and the displaced fluid. The yellow arrow represents the flow direction. The green cylinder represents an impermeable solid material, while the area formed between the cylinder and the cylinder is a pore channel, and the fluid flows in the pore channel between the cylinders.

Model Verification
We first performed mesh independence verification and time step independence verification. The time steps selected for grid independence verification were all 1 × 10 −6 s; the grid size selected for time step independence verification was 6 × 10 −5 m. The total calculation time is the time required for the displacement fluid to flow to the porous media The red area is the displacing fluid, which is injected at a certain flow rate from below. The blue area is the displaced fluid, flowing out from the upper outlet. The outlet gauge pressure is 0 Pa. The horizontal line between the red area and the blue area represents the interface between the displacing fluid and the displaced fluid. The yellow arrow represents the flow direction. The green cylinder represents an impermeable solid material, while the area formed between the cylinder and the cylinder is a pore channel, and the fluid flows in the pore channel between the cylinders.

Model Verification
We first performed mesh independence verification and time step independence verification. The time steps selected for grid independence verification were all 1 × 10 −6 s; the grid size selected for time step independence verification was 6 × 10 −5 m. The total calculation time is the time required for the displacement fluid to flow to the porous media outlet. Figure 4 exhibits the time evolution of the saturation of displacing fluid and the interface length between the displacing fluid and the displaced fluid under different mesh sizes. The results show that the saturation and interface length tend to stabilize when the grid size is 6 × 10 −5 m. Figure 5 presents the time evolution of the saturation of the displacing fluid and the interface length between the displacing fluid and the displaced fluid under different time steps. The results show that the saturation and interface length tend to be stable when the time step is 1 × 10 −6 s. Thus, a grid size of 6 × 10 −5 m and a time step of 1 × 10 −6 s were employed in the subsequent numerical simulations.     We selected the experimental results of Rabbani et al. [11] for experimental verification. Their experiment involves filling silicone oil into a porous medium arranged by cylinders. Then, water is injected into the porous medium to see how water displaces silicone oil in the porous medium. Among them, two different porous structures are used. One is the porous medium with uniform pores. The other is porous media with a pore gradient whose size increases from inlet to outlet. We selected four groups of experimental results to verify the current numerical model, including the fluid displacement processes in the porous mediums with and without pore gradient under the capillary number of 7.5 × 10 −6 We selected the experimental results of Rabbani et al. [11] for experimental verification. Their experiment involves filling silicone oil into a porous medium arranged by cylinders. Then, water is injected into the porous medium to see how water displaces silicone oil in the porous medium. Among them, two different porous structures are used. One is the porous medium with uniform pores. The other is porous media with a pore gradient whose size increases from inlet to outlet. We selected four groups of experimental results to verify the current numerical model, including the fluid displacement processes in the porous mediums with and without pore gradient under the capillary number of 7.5 × 10 −6 and 1.4 × 10 −5 . In order to save computing time and resources, we simulated a part of the experimental model. The fluid parameters and the amplitude of the pore gradient decline are consistent with the parameters in the experiment. The comparisons between simulation results and experimental results are shown in Figures 6 and 7. The uniform medium shown in Figure 6 was composed of cylinders with the same diameter of 1 mm. The nonuniform medium in Figure 6 consisted of cylinders with the diameter decreasing row by row from the top to the bottom. The diameter difference between the adjacent row of cylinders was 0.0135 mm and the diameters of the cylinder in the top and bottom rows were 1 mm and 0.81 mm, respectively.     [11]. results of the current investigation. It can be seen from the comparison between the simulation results an mental results [8] that, when the capillary number is 7.5 × 10 −6 in the nonunif (the left figure of A and B in Figure 6), it can be simulated that the flow of dis in porous media is stable. In the other three cases, it can also be simulated placing fluid has an unstable flow. In order to compare the simulation res It can be seen from the comparison between the simulation results and the experimental results [8] that, when the capillary number is 7.5 × 10 −6 in the nonuniform medium (the left figure of A and B in Figure 6), it can be simulated that the flow of displacing fluid in porous media is stable. In the other three cases, it can also be simulated that the displacing fluid has an unstable flow. In order to compare the simulation results with the experimental results more clearly, we also drew a diagram comparing the saturation under various conditions, as shown in Figure 6. It can be seen from the figure that the simulation results are similar to the experimental results, which shows that the lattice Boltzmann method can simulate fluid flow in porous media.

Effect of Wettability
When the non-wetting fluid displaces the wetting fluid in a porous medium, it will flow forward only when the driving pressure of the displacement flow is greater than the capillary pressure threshold~σcos θ/r, where r is the pore radius. While the capillary pressure threshold can be adjusted by regulating the wettability to further regulate the fluid displacement process. It can be seen from the above formula that, with the increase of wetting angle, the capillary pressure threshold increases and then the flow resistance increases. Based on this, a porous medium with the wettability gradients will hopefully be designed to improve the flow stability, as shown in Figure 8. When the fluid flows from a row of cylinders with a wetting angle of θ to the next row of cylinders with a larger wetting angle of θ', the flow resistance increases and the fluid will preferentially flow transversely and then forward.
wetting angle, the capillary pressure threshold increases and then the flow resistance increases. Based on this, a porous medium with the wettability gradients will hopefully be designed to improve the flow stability, as shown in Figure 8. When the fluid flows from a row of cylinders with a wetting angle of θ to the next row of cylinders with a larger wetting angle of θ', the flow resistance increases and the fluid will preferentially flow transversely and then forward. In order to verify the influence of wettability on the fluid displacement process in the porous medium, the variation curves of average flow rate and outflow time at different wetting angles were investigated, as shown in Figure 9. With the increase in the wetting angle, the average flow rate of the displacement fluid decreases and the time required for the displacement fluid to flow out of the porous media increases. Thus, the flow resistance of the displacement fluid in porous media increases with the increment in wetting angle, which is consistent with the above theoretical analysis results. In order to verify the influence of wettability on the fluid displacement process in the porous medium, the variation curves of average flow rate and outflow time at different wetting angles were investigated, as shown in Figure 9. With the increase in the wetting angle, the average flow rate of the displacement fluid decreases and the time required for the displacement fluid to flow out of the porous media increases. Thus, the flow resistance of the displacement fluid in porous media increases with the increment in wetting angle, which is consistent with the above theoretical analysis results.  Figure 10 shows the simulation results for different capillary numbers and wetting gradients. It can be seen that the flow pattern of the displacing fluid is unstable with or without wetting gradient when the capillary number is large. In the case of wetting gradient, the stability of the displacing fluid flow front is improved when the capillary number decreases. Without wetting gradient, the fluid still has more bifurcation. When the capillary number decreases to 7.5 × 10 −6 , the result with the wetting gradient is the best. The resistance gradient due to the wetting gradient can be fully displayed. In the absence of a wetting gradient, the flow of the displacing fluid is unstable even when the capillary number is reduced. It can be clearly seen from the simulation results in Figure 9 that the wetting gradient plays a significant role in improving the stability of the flow front. Furthermore, the smaller the capillary number (the lower the injection speed), the better the fluid filling, and the more stable the fluid front is.  Figure 10 shows the simulation results for different capillary numbers and wetting gradients. It can be seen that the flow pattern of the displacing fluid is unstable with or without wetting gradient when the capillary number is large. In the case of wetting gradient, the stability of the displacing fluid flow front is improved when the capillary number decreases. Without wetting gradient, the fluid still has more bifurcation. When the capillary number decreases to 7.5 × 10 −6 , the result with the wetting gradient is the best. The resistance gradient due to the wetting gradient can be fully displayed. In the absence of a wetting gradient, the flow of the displacing fluid is unstable even when the capillary number is reduced. It can be clearly seen from the simulation results in Figure 9 that the wetting gradient plays a significant role in improving the stability of the flow front. Furthermore, the smaller the capillary number (the lower the injection speed), the better the fluid filling, and the more stable the fluid front is. To quantitatively analyze the flow results shown in Figure 11, three parameters were selected for the quantitative analysis. The three quantitative parameters were the fractal dimension ( Figure 11A), the displacement efficiency ( Figure 11B), and the interface length ( Figure 11C). It can be seen from the results of Figure 11A that, on the whole, the fractal dimension with the wetting gradient is larger than that without the wetting gradient. This indicates the effect of the wetting gradient on improving the flow stability of displacing fluids; with the wetting gradient, the capillary number is smaller and the fractal dimension is larger. The smaller the capillary number, the better the stability of the fluid front. This also indicates that both the wetting gradient and the capillary number affect the stability of the flow front. From the comparison of the displacement efficiency results of Figure  11B, it can be seen that, on the whole, the displacement efficiency with the wetting gradient is higher than that without the wetting gradient, which indicates that the existence of the wetting gradient improves the filling degree of the fluid. As the capillary number decreases, the displacement efficiency increases, and the filling degree improves. From the comparison of the results in Figure 11C, it can be seen that the interface length of porous media with the wetting gradient is much shorter than that without the wetting gradient, which also indicates that the bending degree of the fluid front interface is reduced. In addition, the smaller the capillary number, the shorter the interfacial length of the fluid front, the smaller the bending degree of the fluid front interface and the better the stability of the fluid front.

Effect of Capillary Number
From the results of these quantitative analyses, when the capillary number increases to 1.8 × 10 −4 under the condition of the wetting gradient, the overall flow stability is reduced and the filling effect is poor. This is because when the capillary number is 1.8 × 10 −4 , the corresponding inlet flow rate is relatively large, and the displacing fluid is subjected to a large force in the vertical direction, which exceeds the threshold range of the resistance gradient. The flow front of the displacing fluid is prone to local preferential flow and, when there is a wetting gradient, due to the effect of the resistance gradient, part of the To quantitatively analyze the flow results shown in Figure 11, three parameters were selected for the quantitative analysis. The three quantitative parameters were the fractal dimension ( Figure 11A), the displacement efficiency ( Figure 11B), and the interface length ( Figure 11C). It can be seen from the results of Figure 11A that, on the whole, the fractal dimension with the wetting gradient is larger than that without the wetting gradient. This indicates the effect of the wetting gradient on improving the flow stability of displacing fluids; with the wetting gradient, the capillary number is smaller and the fractal dimension is larger. The smaller the capillary number, the better the stability of the fluid front. This also indicates that both the wetting gradient and the capillary number affect the stability of the flow front. From the comparison of the displacement efficiency results of Figure 11B, it can be seen that, on the whole, the displacement efficiency with the wetting gradient is higher than that without the wetting gradient, which indicates that the existence of the wetting gradient improves the filling degree of the fluid. As the capillary number decreases, the displacement efficiency increases, and the filling degree improves. From the comparison of the results in Figure 11C, it can be seen that the interface length of porous media with the wetting gradient is much shorter than that without the wetting gradient, which also indicates that the bending degree of the fluid front interface is reduced. In addition, the smaller the capillary number, the shorter the interfacial length of the fluid front, the smaller the bending degree of the fluid front interface and the better the stability of the fluid front.
width [33]. Therefore, we used this method to distinguish whether the flow of fluid in porous media is stable or unstable. It can be seen from the phase diagram that, in the case of the wetting gradient, the displacement process of fluid in porous media is mostly stable except for in the case of a large capillary number. The results of the fluid displacement process are unstable without the wetting gradient, which indicates that the wetting gradient plays a significant role in improving the flow stability of fluids in porous media.  From the results of these quantitative analyses, when the capillary number increases to 1.8 × 10 −4 under the condition of the wetting gradient, the overall flow stability is reduced and the filling effect is poor. This is because when the capillary number is 1.8 × 10 −4 , the corresponding inlet flow rate is relatively large, and the displacing fluid is subjected to a large force in the vertical direction, which exceeds the threshold range of the resistance gradient. The flow front of the displacing fluid is prone to local preferential flow and, when there is a wetting gradient, due to the effect of the resistance gradient, part of the displacement fluid flows laterally, making the bifurcation phenomenon of the displacement fluid front more serious. However, in the case of the wetting gradient, when the capillary number is less than 3.5 × 10 −5 , the overall filling effect is better and the stability of the flow is improved with the decrease of the capillary number. This is because the force of the displacement fluid in the vertical direction is less than the threshold range of the resistance gradient, which means the resistance gradient can inhibit the occurrence of local preferential flow.
We also plotted the phase diagram of the wetting gradient and capillary number, as shown in Figure 12. The identification of flow patterns in porous media can be considered comprehensively by way of parameters such as saturation, fractal dimension and finger width [33]. Therefore, we used this method to distinguish whether the flow of fluid in porous media is stable or unstable. It can be seen from the phase diagram that, in the case of the wetting gradient, the displacement process of fluid in porous media is mostly stable except for in the case of a large capillary number. The results of the fluid displacement process are unstable without the wetting gradient, which indicates that the wetting gradient plays a significant role in improving the flow stability of fluids in porous media. In order to better show the effect of wettability gradient on fluid flow in porous media, we selected a group of typical simulation results for comparison at different times. Figure 13 shows the comparison of saturation and interface length at different times with and without wetting gradient when the capillary number was 1.8 × 10 −5 . The time is normalized by the time required for the displacement fluid to flow to the porous media outlet, to characterize the different stages of the entire displacement process. The slope of the curve in Figure 13A shows that the filling speed of the fluid is faster when there is a wetting gradient, and the saturation result at any normalization time shows that the filling effect is better when there is a wetting gradient than when there is no wetting gradient. Figure 13B shows that the slope increases significantly in the case of no wetting gradient after the normalization time is 0.2, which indicates that there is unstable flow at the front of the fluid at this moment, especially the bifurcated finger-like flow, which makes the interface length of the front of the fluid significantly longer. In the case of the wetting gradient, the slope of the interface length is relatively flat, which means that the flow is relatively stable and there is no bifurcation phenomenon. So, the change of the interface length at the front of the fluid is not very large, which also means that porous media with a wetting gradient can better improve the stability of the flow front.  In order to better show the effect of wettability gradient on fluid flow in porous media, we selected a group of typical simulation results for comparison at different times. Figure 13 shows the comparison of saturation and interface length at different times with and without wetting gradient when the capillary number was 1.8 × 10 −5 . The time is normalized by the time required for the displacement fluid to flow to the porous media outlet, to characterize the different stages of the entire displacement process. The slope of the curve in Figure 13A shows that the filling speed of the fluid is faster when there is a wetting gradient, and the saturation result at any normalization time shows that the filling effect is better when there is a wetting gradient than when there is no wetting gradient. Figure 13B shows that the slope increases significantly in the case of no wetting gradient after the normalization time is 0.2, which indicates that there is unstable flow at the front of the fluid at this moment, especially the bifurcated finger-like flow, which makes the interface length of the front of the fluid significantly longer. In the case of the wetting gradient, the slope of the interface length is relatively flat, which means that the flow is relatively stable and there is no bifurcation phenomenon. So, the change of the interface length at the front of the fluid is not very large, which also means that porous media with a wetting gradient can better improve the stability of the flow front. In order to better show the effect of wettability gradient on fluid flow in porous media, we selected a group of typical simulation results for comparison at different times. Figure 13 shows the comparison of saturation and interface length at different times with and without wetting gradient when the capillary number was 1.8 × 10 −5 . The time is normalized by the time required for the displacement fluid to flow to the porous media outlet, to characterize the different stages of the entire displacement process. The slope of the curve in Figure 13A shows that the filling speed of the fluid is faster when there is a wetting gradient, and the saturation result at any normalization time shows that the filling effect is better when there is a wetting gradient than when there is no wetting gradient. Figure 13B shows that the slope increases significantly in the case of no wetting gradient after the normalization time is 0.2, which indicates that there is unstable flow at the front of the fluid at this moment, especially the bifurcated finger-like flow, which makes the interface length of the front of the fluid significantly longer. In the case of the wetting gradient, the slope of the interface length is relatively flat, which means that the flow is relatively stable and there is no bifurcation phenomenon. So, the change of the interface length at the front of the fluid is not very large, which also means that porous media with a wetting gradient can better improve the stability of the flow front.

Effect of Viscosity Ratio
In order to clarify the influence of wetting gradient on the fluid displacement process with a wider range of parameters, we compared and studied the fluid displacement process of different fluid viscosities with and without a wetting gradient. We chose the case with a capillary number of 3.5 × 10 −5 for our research. This is because under this capillary number, although there is a wetting gradient that plays a certain role in the stable flow of the fluid, there are still bifurcation phenomena at the front of the fluid, which affect the filling effect. Therefore, we tried to study how to improve the flow stability of the front of the fluid by changing the viscosity ratio. Figure 14A shows the comparison results of saturation with or without wetting gradient when the viscosity ratio is 0.0025, 0.005, 0.01, 0.02 and 0.04, respectively. We kept the viscosity of the displaced fluid unchanged at 1 mPa s, and changed the viscosity of the displaced fluid. The viscosities of the displaced fluid were set at 400 mPa s, 200 mPa s, 100 mPa s, 50 mPa s and 25 mPa s. From the comparison results in the figure, it can be seen that the change of viscosity ratio does have a great impact on the flow filling effect. On the whole, the saturation value under the condition of the wetting gradient is higher than that under the condition of no wetting gradient, which indicates that the filling effect under the condition of the wetting gradient is better. In addition, under the condition of the wetting gradient, the greater the viscosity ratio, the better the stability of the fluid front. Figure 14B shows the comparison results of interface length under different viscosity ratios. On the whole, the interface length under the condition of the wetting gradient is smaller than that under the condition of no wetting gradient, indicating that the flow stability under the condition of the wetting gradient is higher. The reason for this result is that when the viscosity ratio increases, the viscosity of the displaced fluid decreases, which makes the flow resistance of the injected fluid relatively lower; for the second half of the fluid flowing into the porous medium especially, the flow resistance decreases to a reasonable size, so that the effect of the resistance gradient can be reflected-the effect of the resistance gradient is to prevent the partial preferential breakthrough when the fluid flows in the porous medium. Therefore, due to the effect of the resistance gradient, the bifurcation phenomenon is reduced. In the case of the wetting gradient, the effect of the viscosity ratio on the wetting gradient is more obvious, so the flow of fluid in the porous media is more stable. In the case of no wetting gradient, the flow resistance of the corresponding injected fluid decreases due to the increase of viscosity ratio. At this time, there is no effect of resistance gradient, and the flow front is more prone to local breakthrough, leading to local preferential flow and affecting the stability of the displacement process.
Materials 2023, 16, x FOR PEER REVIEW 15 of 20 is 1.8 × 10 −5 and the viscosity ratio is 0.005. The solid and hollow symbols indicate the porous media with and without wetting gradient, respectively.

Effect of Viscosity Ratio
In order to clarify the influence of wetting gradient on the fluid displacement process with a wider range of parameters, we compared and studied the fluid displacement process of different fluid viscosities with and without a wetting gradient. We chose the case with a capillary number of 3.5 × 10 −5 for our research. This is because under this capillary number, although there is a wetting gradient that plays a certain role in the stable flow of the fluid, there are still bifurcation phenomena at the front of the fluid, which affect the filling effect. Therefore, we tried to study how to improve the flow stability of the front of the fluid by changing the viscosity ratio. Figure 14A shows the comparison results of saturation with or without wetting gradient when the viscosity ratio is 0.0025, 0.005, 0.01, 0.02 and 0.04, respectively. We kept the viscosity of the displaced fluid unchanged at 1 mPa s, and changed the viscosity of the displaced fluid. The viscosities of the displaced fluid were set at 400 mPa s, 200 mPa s, 100 mPa s, 50 mPa s and 25 mPa s. From the comparison results in the figure, it can be seen that the change of viscosity ratio does have a great impact on the flow filling effect. On the whole, the saturation value under the condition of the wetting gradient is higher than that under the condition of no wetting gradient, which indicates that the filling effect under the condition of the wetting gradient is better. In addition, under the condition of the wetting gradient, the greater the viscosity ratio, the better the stability of the fluid front. Figure 14B shows the comparison results of interface length under different viscosity ratios. On the whole, the interface length under the condition of the wetting gradient is smaller than that under the condition of no wetting gradient, indicating that the flow stability under the condition of the wetting gradient is higher. The reason for this result is that when the viscosity ratio increases, the viscosity of the displaced fluid decreases, which makes the flow resistance of the injected fluid relatively lower; for the second half of the fluid flowing into the porous medium especially, the flow resistance decreases to a reasonable size, so that the effect of the resistance gradient can be reflected-the effect of the resistance gradient is to prevent the partial preferential breakthrough when the fluid flows in the porous medium. Therefore, due to the effect of the resistance gradient, the bifurcation phenomenon is reduced. In the case of the wetting gradient, the effect of the viscosity ratio on the wetting gradient is more obvious, so the flow of fluid in the porous media is more stable. In the case of no wetting gradient, the flow resistance of the corresponding injected fluid decreases due to the increase of viscosity ratio. At this time, there is no effect of resistance gradient, and the flow front is more prone to local breakthrough, leading to local preferential flow and affecting the stability of the displacement process. We selected two groups with a viscosity ratio M of 0.01 and 0.02 for comparative analysis at different times, as shown in Figure 15. On the whole, the saturation with a wetting gradient is higher than that without a wetting gradient, indicating that the filling effect is better with a wetting gradient. In terms of interface length, the interface length with the wetting gradient is smaller than that without the wetting gradient, which indicates that the stability of the fluid front is better with a wetting gradient; after the normalization time is 0.2, the interface length without the wetting gradient increases faster, which indicates that the unstable flow phenomenon begins at this time, and the finger-like flow becomes more and more obvious in the later stage, making the interface length increase faster. At the same time, the saturation of the fluid increases slowly. On the whole, the filling condition with the wetting gradient is better than that without a wetting gradient, and the filling effect with a high viscosity ratio is better than that with a low viscosity ratio. Figure 14. The effect of viscosity ratio on the displacement efficiency (A) and the interface len (B). In both cases, the capillary number is 3.5 × 10 −5 and the porosity is 0.48. The solid and holl symbols indicate the porous media with and without wetting gradient, respectively.
We selected two groups with a viscosity ratio M of 0.01 and 0.02 for comparat analysis at different times, as shown in Figure 15. On the whole, the saturation with wetting gradient is higher than that without a wetting gradient, indicating that the filli effect is better with a wetting gradient. In terms of interface length, the interface leng with the wetting gradient is smaller than that without the wetting gradient, which in cates that the stability of the fluid front is better with a wetting gradient; after the norm ization time is 0.2, the interface length without the wetting gradient increases faster, wh indicates that the unstable flow phenomenon begins at this time, and the finger-like fl becomes more and more obvious in the later stage, making the interface length increa faster. At the same time, the saturation of the fluid increases slowly. On the whole, t filling condition with the wetting gradient is better than that without a wetting gradie and the filling effect with a high viscosity ratio is better than that with a low viscosity rat

Effect of Porosity
The size of the porosity also has a great impact on the flow of fluid in porous med The size of porosity especially affects the function of the wetting gradient. Figure 16 sho the saturation and interface length of different porosities at the normalized time. It can seen that, in the case of a wetting gradient, when the porosity is large, the slope of fluid saturation is large, indicating that the flow front is relatively stable. This is becau in the case of the wetting gradient, with the increase of porosity, the upper and low spacings between cylinders become wider, and the contact area between fluid and poro media decreases. This can reduce the resistance gradient to an appropriate size, and th make the resistance gradient play a maximum role, which inhibits the occurrence of lo breakthrough. On the whole, the filling condition with the wetting gradient is better th that without a wetting gradient, and the larger the porosity, the better the filling effect.

Effect of Porosity
The size of the porosity also has a great impact on the flow of fluid in porous media., The size of porosity especially affects the function of the wetting gradient. Figure 16 shows the saturation and interface length of different porosities at the normalized time. It can be seen that, in the case of a wetting gradient, when the porosity is large, the slope of the fluid saturation is large, indicating that the flow front is relatively stable. This is because in the case of the wetting gradient, with the increase of porosity, the upper and lower spacings between cylinders become wider, and the contact area between fluid and porous media decreases. This can reduce the resistance gradient to an appropriate size, and then make the resistance gradient play a maximum role, which inhibits the occurrence of local breakthrough. On the whole, the filling condition with the wetting gradient is better than that without a wetting gradient, and the larger the porosity, the better the filling effect.

Effect of Nonuniformity
We randomly adjusted the diameter of some cylinders to obtain the nonuniform medium. In order to maintain the same porosity as the uniform medium, we adopted the method of expanding and shrinking the same number of cylinder diameters, respectively, to make the porosity of the nonuniform medium consistent with that of the uniform medium, as shown in Figure 17. In addition, we chose to compare the uniform media with a capillary number of 3.5 × 10 −5 , a viscosity ratio of 0.02 and a porosity of 0.48. From the comparison of saturation and interface length in Figure 18, it can be seen that under the condition of the wetting gradient, the flow stability of the displacement fluid in the nonuniform medium is indeed not as good as that in the uniform medium. This is because the nonuniformity of pores in the inhomogeneous medium causes the displacement fluid to flow preferentially in local areas, which makes the flow stability of the displacement fluid worse. In addition, under the condition of no wetting gradient, because there is no resistance gradient, this local preferential flow phenomenon is more likely to occur, resulting in the flow stability of the displacement fluid in the nonuniform medium being worse. However, from the comparison diagram, even in the case of nonuniform media, the filling effect of the displacement fluid with the wetting gradient is still better than that without a wetting gradient, which also shows the role of the wetting gradient in improving the stability of the displacement process.

Effect of Nonuniformity
We randomly adjusted the diameter of some cylinders to obtain the nonuniform medium. In order to maintain the same porosity as the uniform medium, we adopted the method of expanding and shrinking the same number of cylinder diameters, respectively, to make the porosity of the nonuniform medium consistent with that of the uniform medium, as shown in Figure 17. In addition, we chose to compare the uniform media with a capillary number of 3.5 × 10 −5 , a viscosity ratio of 0.02 and a porosity of 0.48. From the comparison of saturation and interface length in Figure 18, it can be seen that under the condition of the wetting gradient, the flow stability of the displacement fluid in the nonuniform medium is indeed not as good as that in the uniform medium. This is because the nonuniformity of pores in the inhomogeneous medium causes the displacement fluid to flow preferentially in local areas, which makes the flow stability of the displacement fluid worse. In addition, under the condition of no wetting gradient, because there is no resistance gradient, this local preferential flow phenomenon is more likely to occur, resulting in the flow stability of the displacement fluid in the nonuniform medium being worse. However, from the comparison diagram, even in the case of nonuniform media, the filling effect of the displacement fluid with the wetting gradient is still better than that without a wetting gradient, which also shows the role of the wetting gradient in improving the stability of the displacement process.

Effect of Nonuniformity
We randomly adjusted the diameter of some cylinders to obtain the nonuniform medium. In order to maintain the same porosity as the uniform medium, we adopted the method of expanding and shrinking the same number of cylinder diameters, respectively, to make the porosity of the nonuniform medium consistent with that of the uniform medium, as shown in Figure 17. In addition, we chose to compare the uniform media with a capillary number of 3.5 × 10 −5 , a viscosity ratio of 0.02 and a porosity of 0.48. From the comparison of saturation and interface length in Figure 18, it can be seen that under the condition of the wetting gradient, the flow stability of the displacement fluid in the nonuniform medium is indeed not as good as that in the uniform medium. This is because the nonuniformity of pores in the inhomogeneous medium causes the displacement fluid to flow preferentially in local areas, which makes the flow stability of the displacement fluid worse. In addition, under the condition of no wetting gradient, because there is no resistance gradient, this local preferential flow phenomenon is more likely to occur, resulting in the flow stability of the displacement fluid in the nonuniform medium being worse. However, from the comparison diagram, even in the case of nonuniform media, the filling effect of the displacement fluid with the wetting gradient is still better than that without a wetting gradient, which also shows the role of the wetting gradient in improving the stability of the displacement process.

Conclusions
In the present study, the displacement of immiscible fluid in porous media was numerically investigated using the lattice Boltzmann method. Our results demonstrate that the displacement pattern and efficiency can be controlled by the wetting gradient. The flow front stability can be improved by setting a wetting gradient on the porous media, which is confirmed over a wide parameter range of different wetting gradients, capillary numbers, viscosity ratios and porosities. When the capillary number is less than 3.5 × 10 −5 , the viscous fingering can be successfully suppressed and the transition from unstable to stable displacement can be achieved by the wetting gradient. While increasing the wetting gradient has little effect on the fluid displacement process, under the conditions of high viscosity ratio (M > 0.01) and large porosity (Φ > 0.58), the improvement of flow front stability is more significant when adding a wetting gradient. In addition, for the porous media with a wetting gradient, decreasing the capillary number, increasing the viscosity ratio or increasing porosity can improve flow front stability and filling efficiency, which is consistent with that of homogeneous wetting media. The present findings could be helpful for the design of porous products to suppress viscous fingering, which will be of great significance for industrial applications such as composite material preparation and oil exploration.

Conclusions
In the present study, the displacement of immiscible fluid in porous media was numerically investigated using the lattice Boltzmann method. Our results demonstrate that the displacement pattern and efficiency can be controlled by the wetting gradient. The flow front stability can be improved by setting a wetting gradient on the porous media, which is confirmed over a wide parameter range of different wetting gradients, capillary numbers, viscosity ratios and porosities. When the capillary number is less than 3.5 × 10 −5 , the viscous fingering can be successfully suppressed and the transition from unstable to stable displacement can be achieved by the wetting gradient. While increasing the wetting gradient has little effect on the fluid displacement process, under the conditions of high viscosity ratio (M > 0.01) and large porosity (Φ > 0.58), the improvement of flow front stability is more significant when adding a wetting gradient. In addition, for the porous media with a wetting gradient, decreasing the capillary number, increasing the viscosity ratio or increasing porosity can improve flow front stability and filling efficiency, which is consistent with that of homogeneous wetting media. The present findings could be helpful for the design of porous products to suppress viscous fingering, which will be of great significance for industrial applications such as composite material preparation and oil exploration.