## 1 Introduction  
The black oil model is a simplified yet robust approach used in reservoir simulations to represent the phase behavior and flow of oil, water, and gas in subsurface formations. The black oil model plays a crucial role in oil production optimization by simulating reservoir behavior and aiding in decision-making processes [1]. It is widely used in the industry for predicting oil recovery and guiding management strategies [2]. The model's ability to accurately capture natural depletion scenarios makes it a valuable tool for reservoir management. The black oil model is essential for integrating optimization methods with simulation and prediction techniques, enabling the industry to make the optimal decisions regarding oil production [3]. The black oil model primarily simplifies the representation of reservoir hydrocarbons into three components--oil, gas, and water--making it a practical choice for many applications in reservoir engineering and simulation. This model is particularly favored for its computational efficiency and the reduced complexity it offers compared to compositional models, which require the detailed characterization of all hydrocarbon components and their interactions [4, 5, 6]. The black oil model is used to simulate the three-phase flow and predict the behavior of the oil, gas, and water phases under various reservoir conditions [7], and its application ranges from primary and secondary recovery processes to enhanced oil recovery (EOR) operations. Although the black oil model has limitations in accurately capturing the phase behavior of more complex fluid systems, such as gas condensate reservoirs, it remains a cornerstone in reservoir simulations due to its ability to provide conservative calculations for compressible and incompressible multiphase flows [8]. Recent advancements have extended the capabilities of the black oil model to better represent the physical properties of oil-gas mixtures, such as density and viscosity, through the dynamic black oil properties that depend on the fraction of CO\({}_{2}\) in the cell, enhancing its applicability to CO\({}_{2}\) EOR operations [9]. Additionally, efforts to incorporate the effects of large gas-oil capillary pressures and compositional changes in gas injection scenarios have further improved the model's accuracy and robustness [10]. Moreover, the development of sophisticated black oil-based multi-component models for polymer flooding and the integration of black oil data with common equation of state (EOS) models for simulating fluid production in multi-reservoir systems with a common surface network demonstrate the model's adaptability and ongoing relevance in addressing contemporary challenges in reservoir management [11; 12]. The traditional method to solve the black oil model in reservoir simulations is the fully implicit method, which is widely used due to its robustness in handling the complex interactions between different phases and components in the reservoir. This method involves the repeated linearization of large nonlinear systems, which often results in ill-conditioned linear systems that are computationally expensive to solve [5; 13; 14]. The fully implicit approach requires solving a coupled system of equations for pressure and saturation simultaneously, which can be challenging due to the nonlinearity and need for accurate phase behavior modeling, especially during phase transitions [15]. One significant drawback of the fully implicit method is its tendency to produce discontinuities in the discrete system when the phase transitions occur, leading to oscillations or even failure of the Newton iterations used for solving the nonlinear equations [15]. Additionally, the method's computational cost is high, because it necessitates solving large, sparse linear systems repeatedly, which can be particularly burdensome for large-scale reservoir models with high geological heterogeneity [16].