Numerical Simulation Study of Tracking the Displacement Fronts and Enhancing Oil Recovery Based on Ferrofluid Flooding

Water flooding is crucial means to improve oil recovery after primary production. However, the utilization ratio of injected water is often seriously affected by heterogeneities in the reservoir. Identification of the location of the displacement fronts and the associated reservoir heterogeneity is important for the management and improvement of water flooding. In recent years, ferrofluids have generated much interest from the oil industry owning to its unique properties. First, saturation of ferrofluids alters the magnetic permeability of the porous medium, which means that the presence of ferrofluids should produce magnetic anomalies in an externally imposed magnetic field or the local geomagnetic field. Second, with a strong external magnetic field, ferrofluids can be guided into regions that were bypassed and with high residual oil saturation. In view of these properties, a potential dual-application of ferrofluid as both a tracer to locate the displacement front and a displacing fluid to improve recovery in a heterogeneous reservoir is examined in this paper. Throughout the injection process, the magnetic field generated by electromagnets and altered by the distribution of ferrofluids was calculated dynamically by applying a finite element method, and a finite volume method was used to solve the multiphase flow. Numerical simulation results indicate that the displacement fronts in reservoirs can indeed be detected, through which the major features of reservoir heterogeneity can be inferred. After the locations of the displacement fronts and reservoir heterogeneities are identified, strong magnetic fields were applied to direct ferrofluids into poorly swept regions and the efficiency of the flooding was significantly improved.


INTRODUCTION
Ferrofluids are stable colloids composed of small nano-scale solid, magnetic particles coated with a molecular layer of dispersant and suspended in a liquid carrier (Rosensweig, 1997). Thermal agitation keeps the particles suspended, and the coatings prevent the particles from coagulation. Its behavior and property can be controlled and changed by the magnetic field (Odenbach, 2008). Therefore, ferrofluid has many industrial applications, such as dynamic sealing, inertial and viscous damper, magnetic drug targeting, liquid microrobots, and so on (Raj and Moskowitz, 1990;Hou et al., 1999;Scherer and Figueiredo, 2005;Fan et al., 2020).
In the past years, ferrofluid flow in porous media has been the subject of various experimental and numerical studies. Borglin, Moridis, and Oldenburg et al. investigated the flow behavior of ferrofluid in porous media, the results showed that the effect of an external magnetic field on the flow of a ferrofluid in porous media is equivalent to the imposition of a magnetic body force, which means the flow of ferrofluids can be guided Oldenburg and Moridis, 1998;Borglin et al., 2000;Oldenburg et al., 2000). The above works, however, did not consider the impact of the ferrofluid on the magnetic field. The presence of ferrofluids should produce anomalies in the magnetic field because these fluids alter the magnetic permeability. Such anomalies, in principle, can be detected by magnetic anomaly detectors. Based on this feature, Sengupta (Sengupta, 2012) proposed an innovative approach to determining the fracture length and width by injecting ferrofluids into the fracture. When applying an external magnetic field, the ferrofluids in the fracture shall generate a secondary magnetic field which will have a definite phase difference from the original magnetic field. Once the phase difference is detected and processed by software an accurate representation of fractures can be made. Schmidt (Schmidt and Tour, 2012) proposed methods for magnetic imaging of geological structures, the principle is similar to Sengupta. Rahmani et al. used a ferrofluid slug to track the movement of a flooding front and studied reservoir permeability heterogeneity (Rahmani et al., 2014;Rahmani et al., 2015).
Water flooding, as is well recognized, is an effective approach to maintain reservoir pressure and improve oil recovery. However, heterogeneities or fractures present in China's continental reservoirs can significantly decrease the sweep efficiency, leading to poor utilization ratios and low oil recovery. There is still a great potential to enhance the recovery (Han, 1995;Yu, 2000;Han, 2007;Yang, 2009;Han 2010). Therefore, it is urgent to find new oil displacement technology and method to overcome the problem of low oil displacement efficiency caused by reservoir heterogeneity and complex structure. In this paper, we study the potential of using ferrofluid both as a tracer and controllable displacing fluid during the whole flooding process. First, since the superparamagnetic nanoparticles in the ferrofluid can change the magnetic permeability of the flooded region, we study the magnetic anomalies and indicate the locations of displacement fronts, as shown in Figure 1A. And then as a displacing fluid that, guided by strong magnetic fields, can overcome reservoir heterogeneities and improve sweep efficiency and oil recovery, as shown in Figure 1B.

Magnetic Field Calculation
In a quasi-static magnetic problem, the relations between the magnetic field strength H, the magnetic flux density B, the current density J, and the magnetization M are given by Maxwell equations and the magnetic properties of the material of concern: Where μ 0 is the permeability of vacuum. Introducing the vector potential A, which is defined by: with the additional assumption that ∇ · A 0, considering Eqs 3-5 can be transformed into In order to solve the differential Eq. 5, it is formulated into a boundary value problem by using the energy functional (Silvester and K. Chari, 1970;Kamminga, 1975) where Ω is the area of the region of study. The variation δI caused by a small variation of δA is equal to Since the energy functional is stationary to any small variation in A, the variation δI must be zero, which gives In regions occupied by a ferrofluid, the magnetic solid particles in the ferrofluid become polarized with the imposition of an external magnetic field, and the ferrofluid is said to be "magnetized." In this paper, we considered a water-based ferrofluid Hinano-FFW provided by a nanomaterials company in China. Although the composition of the fluid was not revealed, the magnetization curve of the fluid was provided to us. Ferrofluid magnetization curves are generally approximated by simple two-parameter arctan functions of the form : where the subscript ff stands for ferrofluid. For our particular fluid, α 1×10 4 , β 3.5 × 10 -5 . Considering a permanent magnet, and the magnetization M m can be defined by: Where M m is the residual magnetization of a permanent magnet, α m is the angle between the direction of the magnetization of the permanent magnet and the positive direction of the z axis, β m is the angle between the direction of the magnetization of the permanent magnet and the positive direction of the x axis, as shown in Figure 2.
In this paper, we consider the 2-D problem and the quasistatic magnetic problem solved by the finite element method.
The magnetization M m of a permanent magnet can be simplified: Thus, the magnetization term M in the Eq. 8 consists of two parts: The computational region Ω was discretized by triangular grids. The current density J and the vector potential A only have component in z direction: From the requirement that the energy functional takes a stationary value, a set of n equations are obtained: where n is the number of grid point. Substituting Eqs 8, 9, 11, 12 into Eq. 14, we obtain Eq. 15 as follows: The potential in the k-th triangular element is assumed to be equal to where Δ k is the volume of k-th element, A i is the potential at the element's vertex, a i ,b i ,r i are the coefficients associated with the element's vertex positions x i ,y i : Because the magnetic flux density B ∇ × A, considering Eq. 16, the components of magnetic induction can be obtained: Take the Eqs 16, 18 into Eq. 15, we can obtain a set of nonlinear equations: and the solution was obtained by using the Newton-Raphson iterative method. Once the vector potential A was obtained, one can calculate the magnetic flux density by B ∇ × A. Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 759862

Simulation of Ferrofluid Flowing in Porous Media
The magnetization of bulk ferrofluid increases with increasing external magnetic field strength, following Eq. 9. In the case of immiscible two-phase flows in porous media, and additional assumption that magnetization increases linearly with ferrofluid's saturation was given by: When an external magnetic field is applied, the ferrofluid receives a magnetic body force F m μ 0 M∇H. This additional magnetic force term is added to the multiphase Darcy's Equation : where k (k x , k y , k z ) is the permeability tensor, k r is the relative permeability, μ is fluid viscosity, p is the fluid pressure, ρ is the fluid density, and D denotes the depth which is positive on the upward side, g is the gravitational acceleration.
For simplicity, we considered the flow as isothermal and incompressible. Such simplifications are common in the analysis of water flooding. From the law of mass conservation, we know that a fluid in a control volume should meet: (22) where the subscript β (o, w, ff ) distinguishes the fluid, n is the outward normal unit vector of region boundary zV, q mβ is the source term that represents mass change per unit time and unit volume, ϕ is the porosity. Using a finite volume method to discretize the mass conservation equation, and the discrete nonlinear equation was also solved by using the Newton-Raphson iterative method.

EXAMPLES AND DISCUSSIONS
Tracking the Displacement Fonts of Ferrofluid Flooding The first example that we considered is a horizontal porous media model as shown in Figure 3A, the left and right boundary is  injector and producer, respectively. We assumed the top and bottom of porous media are closed boundaries. This model can be used to approximate flooding with a line drive pattern. The porosity is 0.2 and the permeability is 100 mD(10 −3 μm 2 ). A tiny permanent magnet was placed on the left side of porous media to generate a background magnetic field, the dimension is 2 cm × 4 cm and the residual flux density of the magnet is 1.19T, the background magnetic field is shown in Figure 3B.
In the beginning, the model was saturated with oil, and the oil viscosity equals 5 mPa·s and the density equal 850 kg/m 3 . Then we used ferrofluid to flooding this porous media, viscosity and density of ferrofluid are 3 mPa·s and 1,187 kg/m 3 , respectively. The relative permeability of the ferrofluid phase k rff S 2 ff . The relative permeability of the oil phase k ro (1 − S ff ) 2 . Both injection and production rates were set to 0.01 V p /min, where V p is the total pore volume. Since a week magnetic field generated by a tiny magnet, we assumed the magnetic force generated on the ferrofluid was neglected in this case.
The magnetic anomaly is defined as: where |B| current is the norm magnetic flux density during ferrofluid flooding, and |B| background is the normal of the magnetic flux density provided by the magnet before the ferrofluid flooding starts.
As Figure 4 shown, the distribution of magnetic anomalies is basically consistent with the ferrofluid saturation distribution. Near the injector is a high magnetic anomalous area because of the highest ferrofluid saturation, however, the decrease of ferrofluid saturation from the injector to the displacement font leads to a magnetic anomaly decrease. Figure 5 shows the ferrofluid saturation and magnetic anomaly curves from the left side to the right side of the porous media at several flooding periods. Note that the inflection point of each anomaly curve basically corresponds to the displacement fonts. For example, point 1 and point 2 of the anomaly curve at 30 min in Figure 5A, correspond to point 1 and point 2 of the saturation curve at 30 min in Figure 5B. It is worth noting that the area between these two points is exactly the displacement front. Thus, this remarkable fact indicates that one can track the ferrofluid fonts based on the magnetic anomaly information at every displacement stage.

Evaluation of the Major Features of Reservoir Heterogeneity
It can be seen from the above section, that the ferrofluid displacement front can be tracked by capturing magnetic anomalies. Furthermore, the flow path of the ferrofluid can be  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 759862 5 obtained from the displacement front position, thus, the major geological characteristics of the reservoir can be inferred.
Consider a five-spot heterogeneous reservoir as Figure 6A shows, the reservoir size is 100 m × 100 m. The heterogeneous reservoir is divided into four regions with an injection well in the center and a production well at each of the four corners. The injection and production rates were set to 0.01 V p /day. Assuming the reservoir has a porosity of 0.2, the permeability in region 1 is k 1 200 mD, and the permeability in regions 2, 3 and 4 is k 2 k 3 k 4 100 mD. In particular, the second region contains an elliptical low-permeability zone while the fourth region contains an elliptical high-permeability zone. The permeability of the two elliptical zones is 10 and 500 mD, respectively. In the beginning, the model was saturated with oil, and the fluid properties are the same as mentioned above. In order to provide a background magnetic field, an electromagnet EM was placed. The size of EM is 10 m × 0.4 m and the current intensity of EM is 10 5 A/m 2 . The background magnetic field is shown in Figure 6B.   As shown in Figure 7, due to the relatively high permeability of region 1 and region 4, the injected ferrofluid propagates faster in these areas, resulting in larger magnetic anomalies. On the contrary, the magnetic anomalies in region 2 and region 3 are relatively small. In particular, the magnetic anomaly in region 4 presents a fingering phenomenon. This is because there is a high permeability elliptical zone in region 4, which makes the injected ferrofluid flow to the production well faster than the other regions. So that the ferrofluid saturation is high at the same time, and resulting in a large magnetic anomaly. However, because the presence of the low permeability elliptical zone in region 2 blocks the flow of ferrofluid to the production well, the ferrofluid in region 2 has a small spread range and low saturation, so the magnetic anomalies are small. In the same way, the lower overall permeability of region 3 is also the reason for the smaller magnetic anomalies. Figure 8, from the magnetic anomaly curves of production wells in different reservoir regions (star mark in Figure 6A), it can be seen that: 1) In the early stage of production, the magnetic anomalies of production well 1, 4 increased and decreased in production well 2, 3. This phenomenon indicates that the injected ferrofluid flows faster in regions 1, 4 than in regions 2, 3, and reflecting indirectly that the permeability of regions 1, 4 is higher than regions 2, 3; 2) Each magnetic anomaly curve has a rapidly rising time point, which corresponds to the moment when the ferrofluid flows into the production well; 3) Because the distribution of magnetic anomalies reflects the distribution of ferrofluid saturation, the magnetic anomaly change information can be used to analyze the flooding process in the reservoir, and so as to recognize the heterogeneous characteristics of the reservoir in general.

Monitoring-Controlling Integrated Ferrofluid Flooding Technology
From the above studies, it can be inferred that the ferrofluid flooding process can be monitored through the magnetic anomaly information, while the major heterogeneity characteristics of the reservoir can be recognized according to the flooding status. In addition, when there is an external magnetic field, the ferrofluid receives a magnetic force, and its flow behavior can be controlled by the magnetic force.
Consider a heterogeneous reservoir as Figure 9A shows. Assuming the reservoir has a porosity of 0.25, and the permeability is 500 mD. In particular, the reservoir contains two elliptical high-permeability zones. The permeability of the two elliptical zones is both 1,500 mD. In order to generate magnetic anomalies during the flooding process, an electromagnet EM1 is placed near the eastern part of the reservoir to provide a background magnetic field. The  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 759862 electromagnet EM1 with a size of 5 m × 0.4 m and a current density of 10 5 A/m 2 . The background magnetic field is shown in Figure 9B. It is worth noting that the background magnetic field is a weak magnetic field, and the generated magnetic force in ferrofluid is very small. So that its influence on the flooding can be ignored. The reservoir is initially saturated with oil, and ferrofluid is used to displace the reservoir. The fluid properties are the same as mentioned above. The injection and production rates were set to 0.025 V p /day. As the high-permeability zone is the dominant channel for flow, the injected ferrofluid mainly flows along the highpermeability zone, so the swept areas are mainly concentrated in the central and southern parts of the reservoir. When the 1.5V p ferrofluid was injected, the oil in the central and southern areas of the reservoir was almost completely displaced, but the oil in the northern area of the reservoir was not displaced, as shown in Figures 10A-D. The distribution of magnetic anomalies in Figures 10A9-D9 shows: 1) At the early stage of flooding, the magnetic anomaly in the central region of the reservoir is relatively large, while the magnetic anomalies in the eastern and northern region of the reservoir are small; 2) At the middle stage of flooding, the magnetic anomaly in the central region of the reservoir is relatively large. These areas gradually increase and spread to the northern region of the reservoir; 3) At the later stage of flooding, the magnetic anomaly in the central and southern regions of the reservoir are generally large, while the magnetic anomaly in the northern region is still low.
From the analysis of the magnetic anomaly changes in the entire flooding process, it can be known that due to the existence of the high permeability zone, the injected ferrofluid is mainly displaced in the central and southern parts of the reservoir, so that the saturation of the ferrofluid is larger, and the magnetic anomaly is larger. Figure 11 shows the magnetic anomaly curves of different monitoring spots in the reservoir, and the positions of the monitoring spots are shown in Figure 9A (star mark). It can be seen from the curves that the ferrofluid mainly flows to the production well and the eastern part of the reservoir because of the existence of the high permeability zone. Therefore, the magnetic anomaly curve of S1 and S2 which is located in the production well and the eastern part of the reservoir gradually increases with the flooding process, while the magnetic anomaly curve of S3 which is located in the northern part of the reservoir basically remains unchanged.
The magnetic anomaly curves of the three monitoring spots show that the injected ferrofluid is more inclined to flood the central and southern areas of the reservoir. Thus, there may be a "dominant flow channel" in this area. The northern part of the  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 759862 8 reservoir has not been affected, and there is still remaining oil in this area, which is the key area for development.
According to the analysis of the above magnetic anomaly curves, it can be seen that the ferrofluid displacement is mainly in the central and eastern parts of the oil reservoir. Therefore, the electromagnet EM2 is deployed near the northern part of the reservoir from the 120th day, and the strong magnetic field generated by EM2 is used to apply a huge magnetic field force on the ferrofluid, leading the ferrofluid to displace to the northern part of the reservoir. The electromagnet EM2 with a size of 10 m × 0.4 m and a current density of 10 8 A/m 2 , and has a running time of 200 days. Figure 12 shows the ferrofluid saturation and magnetic anomaly distribution in the process of ferrofluid flooding after the application of a strong magnetic field. It can be seen that under the control of the external strong magnetic field, the injected ferrofluid can overcome the influence of the reservoir heterogeneity and displace towards the northern part of the reservoir, so that to recover the remaining oil in these areas.
As shown in Figure 13, the magnetic anomaly curve of the monitoring spot S3 in the northern part of the reservoir began to rise rapidly (the dotted red line is without a strong magnetic field) after EM2 was applied. This reflects the rapid increase of ferrofluid saturation, meaning that the ferrofluid is guided to displace to the northern part of the reservoir.
In addition, the water flooding is simulated for the same reservoir. Due to the large difference in oil and water viscosity, the injected water is easier to flow along the high permeability zone to the production well. Therefore, the scope of water flooding in the reservoir is smaller, and there is a large amount of remaining oil in the northern area of the reservoir that has not been displaced, as shown in Figure 14. After 600 days of production, the oil recovery of water flooding is 63.3%, and the oil recovery of ferrofluid flooding under the control of an external magnetic field is 79.5%.

CONCLUSION
In this paper, a novel ferrofluid flooding method was proposed. The ferrofluid plays the dual role of "tracer" and displacement fluid in the flooding process, that is, the flooding process can be monitored and the flooding path can be controlled.
The simulation results show that: 1) By monitoring the magnetic anomaly change information in different areas of the reservoir, the flow direction of the injected ferrofluid can be reflected, so that the heterogeneous characteristics of the reservoir can be understood; 2) The displacement direction can be changed artificially by applying a strong magnetic field to further improve the sweep range; 3) Compared with the traditional water flooding, the monitoring-controlling integrated ferrofluid flooding technology can greatly improve oil recovery.
However, there are still many problems and challenges in applying ferrofluid flooding technology. For example, how electromagnets are placed in the field and the retention loss of ferrofluid during flooding process. In the future, the authors will develop the ferrofluid flooding numerical simulation method for three-dimensional problems, and carry out ferrofluid flooding experimental research to further verify the feasibility of this technology. Nevertheless, the research in this article can still provide new ideas and theoretical guidance for further enhancing oil recovery.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
TH conceived the idea and wrote the manuscript. FS, RW, and XH contributed to the study and gave final approval for the manuscript.  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 759862