2.1 Fast, accurate, and robust phase equilibrium model

In the study of multicomponent and multiphase flow, especially when the complex porous structure and thermodynamic conditions (temperature, pressure) of deep reservoirs are considered, one key step is determining the phase composition and distribution of the fluid mixture that flows, that is, determining whether it is a single gas phase, a single liquid phase, a gas–liquid mixed phase, or a supercritical phase and whether phase transition is underway (Zhang et al., 2024; Zhang, Zhang, et al., 2022). Based on this, other physical and chemical properties of the fluid can be calculated, such as its density, the composition of each component phase, and its surface tension. This process is also known as flash calculation (Zhang et al., 2021). Flash calculation can usually be divided into two categories: NPT (given pressure and temperature) flash and NVT (given molar volume and temperature) flash algorithms. The NPT flash algorithm has a long history of development and simple steps and is now widely used in the energy industry (Babayemi & Eluno, 2021). However, it has been noted that the NPT flash method has two significant drawbacks. First, the results of NPT flash may not provide a unique solution to a cubic equation, which requires human intervention to select physically meaningful results and poses obstacles to large-scale flash calculations (Michelsen, 1982). There may be extreme cases where there are virtual roots or where the roots have no physical meaning, which also hinders the application of this method in engineering research. In addition, another key drawback of the NPT flash algorithm is that in actual multicomponent reservoir fluid flow processes, pressure often cannot be given as a prerequisite, but is a thermodynamic variable that needs to be calculated (Zhang et al., 2023). This also significantly limits the practical application of the NPT flash algorithm in engineering practice. Thus, NVT flash is developed in this section with energy stability properties and adheres to the second law of thermodynamics.

The traditional energy industry is expanding production capacity and upgrading technology with the rapid development of artificial intelligence and big data, and the capacity for computing reservoir simulation often reaches the level of billions of grids (Al-Saadoon et al., 2013). For the common multiphase flow problems in reservoir simulations, phase equilibrium analysis must be carried out simultaneously in these billions of grids to clarify the existence of multiphase flow and provide a physically meaningful initial field. The iterative flash calculation algorithm cannot fulfill the needs of such large-scale and high-precision calculations, and thus more advanced modeling and simulation methods are needed to achieve an acceptable computational cost in engineering.

2.2 FSMA algorithm

To address the pore network extraction problem, we proposed a faster and more accurate FSMA method (Liu et al., 2023). The minimum distance from a void space to a solid space can be calculated using the function in Figure 1. We use this function to assign each point in the void space a dist value, which is used to identify the pore center or critical point. Initially, one point in the void space is randomly selected, and the steepest descent method is used to identify the first pore center. The search space around this pore center (the red point) can be discretized, as in the search region in Figure 1, and the dist value of each discretized point is then determined. Then, along the green circle (the search region), the sequence of discretized points forms a distal curve.

As we introduced in our previous study (Liu et al., 2023), the discontinuous point of the distal curve is the medial axis point, which is named the critical point. Thus, the critical points are determined in the search region, which is shown as three blue points in Figure 1. These three critical points from the pore center represent three different branches of the pore network. By following the different branches, the entire pore network can be determined. As a consequence, the problem becomes critical point identification from one pore center to the next pore center. Inspired by the flashlight, we try to search for the critical point using a fan-shaped search region in front of the previous point.

On the curve of the fan-shaped search region, we use a similar idea. The curve is discretized and the dist distribution is obtained. Accordingly, the critical point is found within the search region. When we repeat the same processes, the identified critical point sequence can represent the medial axis of the void space. After the sequence of critical points is obtained, the dist value sequence is shown. On the dist distribution, we can observe the local maximum and local minimum dist values, which correspond to the pore center and the throat center, respectively.

2.3 Fully conservative implicit pressure explicit saturation (IMPES) algorithm

Darcy-scale modeling and simulation of reservoir porous media have undergone considerable development, and conventional methods have been proposed to describe the flow and diffusion phenomena of gases and fluids in solid porous media. These mathematical models tend to have strong nonlinearity. There are usually three solution schemes for nonlinear mathematical models: fully implicit, semi-implicit, and fully explicit. The fully implicit processing of all items and the simultaneous resolution of all unknown items are therefore often considered an unconditionally stable algorithm. However, implicit methods require significant computational resources for each time step. The IMplicit–EXplicit (IMEX) format implicitly processes linear terms and explicitly calculates other terms (Gardner et al., 2018), and it has better stability than the fully explicit format and can eliminate the nonlinearity of the original equation. Compared with fully implicit algorithms, it has better computational efficiency. Therefore, the IMEX format has been used in the design of numerical simulation algorithms for diffusion layer systems. Operator decoupling is an effective method that simplifies complex time-dependent problems into several simpler ones. Using operator decoupling, iterative operator splitting schemes can be further constructed and developed into iterative methods for solving nonlinear systems. The IMPES scheme has been widely used to solve coupled nonlinear systems of a two-phase flow in porous media (Chen et al., 2004). It is a typical IMEX format that uses a physics-based decoupling method. In the standard IMPES format proposed in the 1960s, the pressure equation was obtained by substituting saturation constraints and Darcy's law in the sum of the two conservation laws of mass for each phase at a continuous level in time and space and by explicitly processing all other variables in the pressure equation to eliminate nonlinearity. Then, the pressure equation was implicitly solved, and as long as the pressure was obtained, the Darcy velocity and two-phase saturation were explicitly updated. Compared with the fully implicit scheme, this method is simpler to set up and more efficient to implement, and it requires less computer memory at each time step.

Because the IMPES scheme has been used for over half a century, its advantages have been recognized and widely accepted by researchers around the world. The IMPES formula is more convenient to implement than the fully implicit scheme. Because of the weak dependence of pressure on saturation and the slow change in pressure over time, the IMPES formula successfully decouples pressure from saturation. Meanwhile, we note some drawbacks of the classical IMPES format: it requires a relatively small time step for numerical stability (as well as to ensure the positive definiteness of saturation). For non-wetting phases, this scheme is not (locally or globally) conserved. Therefore, for the whole fluid mixture, this algorithm is not conservative, locally or globally. In addition, the algorithm may generate a wetting phase saturation greater than 1, which renders its results physically meaningless. In the presence of discontinuous capillary pressure, it cannot reproduce the correct saturation distribution with discontinuity. According to the literature and a theoretical analysis, capillary pressure should play a crucial role in the transport flow of the gas–liquid two-phase flow in deep oil and gas reservoirs; therefore, accurate modeling and description are needed. To address the shortcomings of the classic IMPES algorithm, Hoteit and Firoozabadi proposed an improved IMPES algorithm (Hoteit and Firoozabadi, 2008). With the H–F version of the IMPES algorithm, they hope to accurately manage the differences in capillary pressure in heterogeneous permeable media, which would have a significant impact on the flow development trend of two-phase immiscible flows. The advantages of the H–F IMPES algorithm have been recognized, mainly its ability to reproduce saturation distributions with the expected discontinuities for different capillary pressure functions. However, similar to the standard IMPES, this algorithm is not (locally or globally) conserved for non-wetting phases. Therefore, for the entire two-phase flow fluid, this scheme is not (locally or globally) conserved. More importantly, similar to the standard IMPES, this algorithm may generate a wetting phase saturation greater than 1.

Based on this, the Sun Algorithm can be described as given at time step , we seek the solutions at time step as follows:

Step 1. Seek as a function of ,

3.1 Phase equilibrium analysis

Usually, as the burial depth increases, the percentage of gas reservoirs also increases. Considering that most of the gas components discovered thus far include representative components such as hydrogen, carbon dioxide, and methane, we first designed a natural gas mixture, as shown in Table 1, for phase equilibrium analysis. The proportion of methane, carbon dioxide, and hydrogen is 0.90:0.01:0.09 for Case 1 and that for Case 2 is 0.09:0.90:0.01. Case 1 represents a natural hydrogen reservoir that is typically generated in deep layers (Zgonnik, 2020) and Case 2 represents the carbon dioxide sequestration scenario in deep reservoirs.

The effect of temperature on phase equilibrium conditions was further analyzed by investigating the temperature change process. The number of phases was calculated at different temperatures with different overall mole concentrations, and the results are plotted in Figure 2. In general, the fluid mixture for Case 1 changes from the vapor–liquid two-phase region to the single-vapor-phase region with a temperature change from 150 to 300 K at all sampled mole concentrations. If the overall mole concentration is larger, which indicates that there are more gas molecules at a specific volume with a higher pressure, a higher temperature is required to enter the single-phase region. For example, 150 K is needed to enter the single-phase region at an overall mole fraction of 200 mol/m3, but this number increases to 156 K if a higher mole fraction of 1000 mol/m3 is fed to the domain. The critical temperature reaches 172 K for a very high mole concentration of 2000 mol/m3, which suggests that a relatively larger temperature is needed for the deep reservoir to maintain the gas phase storage reserve for a long time. Similarly, if we want to improve the production rate and transportation flux of natural hydrogen reservoirs, a phase stability test should also be carried out to check whether the temperature needs to be increased to avoid liquefaction during the process.

The number of phases at equilibrium that change with the temperature at different overall mole concentrations for Case 2 is plotted in Figure 3. Compared with Figure 2, a similar change from the vapor–liquid two-phase region to the single-vapor-phase region with temperature can be observed, and a larger critical temperature is always required for a higher overall concentration. In addition, the critical temperature for Case 2 is higher than that for Case 1 at the same overall mole fraction, which is reasonable and can be explained by the higher critical temperature of the dominant carbon dioxide composition. For example, at an overall mole concentration of 200 mol/m3, the operational temperature should reach 208 K to ensure a single vapor phase in the fluid mixture for Case 2, whereas the entering temperature is only 144 K for Case 1. Moreover, the difference in the critical temperature at different overall mole concentrations for Case 2 is larger than that for Case 1. For example, the entering temperature increases from 208 to 240 K with concentrations ranging from 200 to 1000 mol/m3 for Case 2, while it increases from 144 to 156 K for Case 1. In other words, a higher temperature increase should be introduced for Case 2 to ensure the single vapor phase if the storage density and transportation flux are improved.

The prediction performance of the deep learning model is shown in Figure 4. The reliability of the deep learning algorithm is validated because its prediction, which is plotted using symbols, matches the ground-truth data, that is, the iterative flash calculation plotted by lines. The trained deep learning model can accurately capture the phase-splitting phenomena by predicting the compositional mole fraction of the components in the vapor phase for Case 2. In particular, when the pressure increases from 50 to 350 MPa, more carbon dioxide is liquefied, which continuously decreases the mole fraction of carbon dioxide in the vapor phase. Meanwhile, the liquefaction rate of methane decreases; thus, the mole fraction of methane in the vapor phase increases slightly to become even higher than that of carbon dioxide. Furthermore, the fraction of hydrogen increases more considerably, which is also reasonable because hydrogen has the lowest critical temperature and pressure and thus is the most difficult to liquefy among the three components. Using deep learning predictions, the energy industry can rapidly obtain basic knowledge about the thermodynamic conditions of production fluids in deep reservoirs, which could help facilitate the process of decision-making for reservoir monitoring, purification, storage, and transportation. For example, if unwanted liquefaction occurs under certain conditions as predicted by the deep learning model, the engineer may reduce the pressure or improve the temperature condition to keep the reservoir production mixture in the single-vapor-phase region. The validity of the proposed deep learning algorithm has been proven to be reliable in extreme conditions, including under very high pressures such as 350 MPa, which describes the specific scenarios of deep reservoirs. Furthermore, as thermodynamic properties are involved in the network architecture, the components in the mixture may vary with different scenarios and may be represented well using thermodynamic features, which reduces the cost of repeated retraining.

3.2 Porous structure generation

The porous medium model was constructed randomly, as shown in Figure 5. The porous model was sealed by using a solid box, which provided the enclosed pore space. The size of the porous medium was 1000 × 1000, where the numbers of solid phase points and void phase points were 580 734 and 419 266, respectively. Thus, the porous medium matrix is sparse, and it enables the acceleration of the computational process, as the FSMA process only needs the neighbor information from the solid phase to calculate the dist.

In the porous medium, the void space is composed of bulk voids and narrow channels, which are called pore throats. With the pore network modeling (PNM) method, the bulk void is considered to be the pore body, which plays a role in storing reservoir fluid. The pore bodies are connected via the throat channel. As shown in Figure 5, the pore centers were captured accurately by the FSMA method and the pore throats were identified accordingly. In this case, we ignored the blind branches in the dead-end corner to simplify the pore network, which could then accelerate the following PNM calculation. However, apart from the dead-end corner, some pore bodies that are disconnected could also be regarded as blind branches, but this depends on the specific problem, such as the mobility of residual oil in the reservoir (Liu et al., 2022).

To validate the resolution independence of the FSMA method, three cases with different resolutions were examined, as shown in Figure 5. The higher resolution model was constructed using the Kronecker product of the original matrix (1000 × 1000) and a matrix of ones. In different resolution scenarios, the medial axis can always be accurately extracted, and the pore and throat centers are in good agreement with each other. Higher resolution models induce more solid phase points and require greater computational resources. Using the appropriate neighbor-searching method helps to accelerate the extraction process, and we adopt matrix-based neighbor-searching to determine the neighbor solid phase within the search range.

Different cases of porous media were constructed to validate the feasibility of the FSMA method in different situations. In deep underground reservoirs, the rock is highly compressed due to the overburden pressure (Luboń & Tarkowski, 2023; Molíková et al., 2022). As a result, the rock porosity is lower than that in a regular reservoir. Tiny throat channels are usually ignored because they present high resistance to fluid flow, but spontaneous imbibition is also important in relation to small pores or throats. As depicted in Figure 6, small throats and pores can be reached; otherwise, these pores would be considered isolated pores that have no contribution to fluid flow. In all the results, the pore network structure in different porous media can also be accurately described.

3.3 Darcy-scale two-phase flow in porous media

To validate the conservation property of the proposed Sun Algorithm, a test case is designed to investigate the Darcy-scale two-phase flow in porous media. A constant rate of wetting phase injection is posed in the left center of one REV (representative element volume [Zhu et al., 2013]) to replace the non-wetting phase in the oil/gas production domain. The saturation distribution of the wetting phase in the replacement process was calculated using the developed Sun Algorithm, and the results are plotted in Figure 7. As can be seen from the saturation profile, the wetting phase is transported from the left inlet to the right outlet, and the non-wetting phase is replaced and produced from the right outlet. The mass conservation property is shown in Figure 8. The total wetting phase in the domain constantly increases at the beginning, as constant injection is ensured and no wetting phase is output outside the domain. After 0.62 years, the front of the wetting phase reaches the right outlet and is ejected from the domain. Thus, the increase in the total wetting phase curve slows down and the wetting phase output curve begins to increase. The conservation property is proven by the equivalence between the total injection of the wetting phase and the sum total of the wetting phase that remains in and that is output from the domain.

The robustness of the Sun Algorithm can be further validated by a numerical test in extreme cases. In practical deep reservoir development scenarios, the density ratio between the wetting and non-wetting phases may be high, especially when the process of carbon dioxide sequestration in deep saline layers is considered (Rosenbauer and Thomas, 2010). To show that the Sun Algorithm can be widely adapted in complex engineering scenarios, a case with a density ratio of 1:1000 is designed and the inlet position is modified. As shown in Figure 9, the wetting phase still flows from the left inlet to the right outlet, but the flow routine and saturation distribution are significantly different. However, the displacement process in the domain can still be stably simulated using the Sun Algorithm, and the changing saturation profile is physically reasonable and reliable in providing production information for the petroleum industry.