Systematic comparison of advanced models of two- and three-parameter equations to model the imbibition recovery profiles in naturally fractured reservoirs

The imbibition process can be considered one of the most important mechanisms during the production of naturally fractured formations. In this process, the hydrocarbon production in the matrix blocks surrounded by water-filled fractures can be described by a recovery curve. Different equations have been proposed to describe the recovery process during the imbibition. This work presents a detailed analysis of the performance of two- and three-parameter models including Weibull, Probit, Logit-Hill, one-hit-multi-target (1HMT), and all-hit-multi-target (AHMT). These models were tested against different experimental and numerical simulation data in a wide range of rock and fluid properties and matrix dimensions. Particularly, the functionality of these models was examined in early, mid, and late times. It should be highlighted that the three-parameter models (1HMT and AHMT) have not been used previously to describe the imbibition data. The results show that the three-parameter models are more accurate to describe the imbibition recovery trends compared to the two-parameter models. Moreover, the analysis revealed that the AHMT model is better for the early-time data (Error = 0.5), the Logit-Hill model is more accurate for the mid-time data (Error = 0.075), and the Weibull model can best fit the late-time imbibition data (Error = 0.04). Finally, the best model for predicting the recovery factor in fractured reservoirs is model 1HMT because the lowest average RMSE (Root-Mean-Square Error) value of 0.0165 was obtained. The findings of this work can be used to more precisely select the model to curve fit the imbibition data.


Abbreviations
1HMT One-hit-multi-target AHMT All-hit-multi-target IOR Improved oil recovery NFR Naturally fractured reservoir RF Recovery factor RMSE Root-mean-square error

Introduction
Considering the importance of energy supply in human daily life and the prominent role of hydrocarbons and their impact in this vital matter, it is important to know and study the oil and gas reservoirs in detail. Meanwhile, naturally, fractured reservoirs (NFRs) are known as important reservoirs with different behavior compared to conventional (unfractured) reservoirs (Aguilera 1980;van Golf-Racht 1982;Saidi 1987). These heterogeneous reservoirs include two separate systems of matrixes and fractures, which are separated from each other as shown in Fig. 1. Matrix blocks are the original reservoir rock before natural fractures occur. Matrix blocks are known for their unique permeability and porosity, as well as fractures by their permeability and porosity. There is a mass transfer between low permeability and high porosity matrix blocks with high permeability and low porosity fracture network. For this reason, this category of reservoirs is called double permeability and double porosity reservoirs (Firoozabadi 2000;Rose 2001;Qasem et al. 2008).
Determining the characteristics and predicting the production behavior of NFR are known as one of the most challenging issues in the oil industry (Ghaedi et al. 2015). The function of predicting the behavior of NFRs production requires basic knowledge of various mechanisms that occur in this type of reservoir. These mechanisms include water imbibition, molecular diffusion, convection, and gas-oil and water-oil gravity drainage (Ardakany et al. 2014;Abbasi et al. 2016). The main difference in recovery and production performance between NFR and non-NFR is due to the difference in capillary pressure in the matrix and fracture environments (Babadagli 1996;Bentley 2022;Tian et al. 2023).
The displacement of oil by water occurs in the initial recovery by mechanisms resulting from aquifers and active waters that affect the production process from reservoirs; then, in improved oil recovery (IOR) methods (Bassir et al. 2023), displacement of the remaining oil in the reservoir is achieved by waterflooding (Jing et al. 2017). When the reservoir rock is water-wet, the process of displacement of oil by water is called spontaneous imbibition. In the case that the displacement of oil by water is accompanied by external energies, it is called forced imbibition, which causes the creation of an additional pressure gradient. The force resulting from the imbibition process is in interaction with other types of forces such as gravity, viscosity, and capillarity forces (Cheng et al. 2022;Zuo et al. 2023). This mechanism has wide applications in various sciences including petroleum engineering, groundwater engineering, engineering geology, soil physics, and civil engineering (Zhang and Austad 2006;Standnes 2010a;Mirzaei-Paiaman et al. 2011a).
The first model introduced to predict oil recovery behavior in the imbibition process was the Aronofsky model (Aronofsky et al. 1958). This model is based on the following assumptions: (1) There is a perfectly continuous mathematical relationship between the recovery factor and time, in which the recovery factor converges to a fixed number with increasing time.
(2) All properties that affect the output rate of the cores are assumed to be constant during the process and do not change.
The Aronofsky model can be represented as follows: where RF(t) is the oil recovery factor as a function of time, RF max is the ultimate oil recovery factor, t is imbibition time, and a is the parameter that must be fitted. Although the Aronofsky model is very simple, the accuracy of this model in certain circumstances is not sufficient.
Given this problem with Aronofsky, (Shouxiang et al. 1997) proposed a generalized equation by using dimensionless imbibition time: (1) where t D is dimensionless time, k is absolute permeability, φ is fractional porosity, is surface tension between oil and water, µ is the viscosity of water and oil, respectively, and Lc is the characteristic length (Mirzaei-Paiaman et al. 2011b).
where V b is the bulk volume of the matrix block, A i is an area of the ith imbibition surface, and L i is the length from the ith imbibition surface to the no-flow boundary. This equation is useful in all rock and fluid conditions and especially in fitting the imbibition process equations in strongly waterwet samples. The fit parameter was approximately a = 0.05. The second model presented to describe and investigate the recovery factor mechanism is Lambert's W function. (Fries and Dreyer 2008) considered flow through capillary tubes and showed that it is possible to solve the Washburn equation for vertical flow including proportional gravity height (Washburn 1921).
where r is tube radius, ρ is water phase density, θ is contact angle, and g is the acceleration due to gravity. If Eq. (5) is divided into L (length of the capillary tube), b cL the term (e.g., saturation) is known as the capillary rise to gravity head. Then, the equation can be written as follows: Second order derivative Second order logarithmic derivative Second order logarithmic derivative It is very suitable for water-wet samples and the inclusion of gravity force. The Lambert's W function is defined as below: In this domain, the maximum relative error is about 0.1%. e is Euler's number: e = 2.718283 … (Ghaedi and Ahmadpour 2020) performed single and two-parameter models for fractured reservoirs in the imbibition process. They were able to get good accuracy from the models they used. The Dose-response model is a biological model that shows the effect of a drug on a specific microbe (Peleg et al. 1997;Toth et al. 2013;Demidenko et al. 2017). Because the shape of this diagram is quite similar to the recovery factor diagram, we used it in the simulation. Here, the mortality rate corresponds to the recovery factor and toxin concentration corresponds to time.
It is important to model a key process such as the spontaneous imbibition process in NFRs. Thus, it is very important to develop equations to model the imbibition tests for NFRs with high accuracy in the prediction of oil production from the matrix blocks. In past studies (Ghaedi and Ahmadpour 2020) and (Standnes 2010b), simple, one-parameter, and two-parameter models were used. The main focus of the current study is to use two-and three-parameter models, to curve fit the production recovery curves for the imbibition process. For this purpose, statistical parameters such as R 2 and RMSE (Root-Mean-Square Error) are used to compare Second order derivative Second order logarithmic derivative Second order logarithmic derivative  Simplified equation Weibull Logit-Hill the presented models with the real data and laboratory experiments. Furthermore, the performance of these models is examined in the early, middle, and late periods.

Model definition and mathematical analysis
In this study, five different mathematical models were investigated to find the best model for describing the water imbibition process in strongly water-wet reservoirs. Models such as Weibull, Probit, and Logit-Hill were evaluated as two-parameter models, and One-Hit-Multi-Target (1HMT) and All-Hit-Multi-Target (AHMT) models were examined as three-parameter models (Demidenko et al. 2017 Weibull dose-response model (also known as the Gumbel model) is one of the most common and widely used toxicity models. This model uses Poisson distribution to estimate different parameters (Carlborg 1981;Christensen 1984; Christensen and Chen 1985). Table 1 shows the Weibull model's main equation and derivatives.
The Probit model was used by Finney in (1971) to estimate the probability of mortality to the concentration of toxic (Finney 1971). Table 2 shows the Probit model's main equation and derivatives.
The Logit-Hill dose-response model was first presented by Berkson in (1944). This model, which is one of the oldest models used in epidemiology, shows a mathematical description of the relationship between the amount of toxic concentration and the mortality rate. Another application is in the use of toxicology to estimate various parameters (Berkson 1944). Table 3 shows the Logit-Hill model's main equation and derivatives.
The multi-target model is adapted from the cell survival models presented by (Casarett 1968;Yakovlev et al. 1993;Bampoe et al. 2000). The main equation and derivatives of 1HMT and AHMT models have shown in Tables 4 and 5, respectively.  Fischer et al. (2006) For real data, the semi-log plot will result in an S-shaped function that exhibits important mathematical behavior. To have a consistent/stable description of recovery behavior in the reservoir, the proposed models must have the same properties. These properties are listed below: 1. As t tends to 0, the value of the function tends to 0.
Because, at the beginning of the water imbibition process, there is no recovery.
2. As t tends to +∞ , the value of the function reaches a maximum value ( rf max ) . At the end of the process, recovery reaches a maximum value which usually is not 1.
3. At both ends (early and late), the change in the S-shaped function (in a logarithmic scale) will not change significantly. In other words, logarithmic first derivatives of the function tend to be 0 at both ends.
4. In a logarithmic scale, the function is monotonically ascending. In other words, the logarithmic first derivative of the function is always positive.
> 0 5. At a middle point (a positive value of t), the direction of the concavity of the logarithmic scale function changes. In other words, the logarithmic second derivative has one positive root.
As mentioned earlier, among all five models, three models are two parametric models (Weibull, Probit, and Logit-Hill) with parameters a and b, and the rest are three parametric models (1HMT, AHMT) with parameters a, b, and p. As can be seen in Table 6, on the one hand, the two parametric models show constant recovery factors in their inclination points. This leads to errors in the predicted recovery values, especially in the middle of the imbibition process. Moreover, in three parametric models, the value of the recovery factor in the inclination point depends on the value of a third parameter (p). For example, the constant values of 1 2 for two models Probit and Logit-Hill and the constant value of 1 − e −1 for the Weibull model. Also, for three-parameter models, the values of 1 − p p+1 p and p p+1 p for 1HMT and AHMT models, respectively, are dependent on the value of the third parameter, i.e., the p parameter. This leads to better predictions in comparison with two-parameter models. Also, by performing algebra and using foresaid mathematical constraints, we found model parameter constraints for each of the models. These numerical conditions are then used for fitting models into datasets. Table 6 shows these parametric conditions and inclination point values. Table 7 shows the simplified model equations.

Data
The laboratory data presented in the previous studies were used in this work to simulate the imbibition process (Morrow and Xie 2001;Zhou et al. 2002;Fischer et al. 2006;Standnes 2010a;Ghasemi et al. 2020). These data also include recovery factor data in oil and gas reservoirs in strong and weak water-wet rocks. In Table 8, different properties such as permeability, porosity, and oil and water viscosities are shown which cover a wide range of rock and fluid properties. Figures 2, 3, 4, 5 and 6 show the recovery factor graphs related to the works of Morrow and Xie (2001), Ghasemi et al. (2020), Standnes (2010a), Fischer et al. (2006), and Zhou et al. (2002). Finally, it should be noted that the recovery imbibition curves of Standnes (2010a) were obtained for strong water-wet rocks d 2 f d(ln (t)) 2 t r = 0fort r > 0 Fig. 6 Recovery curves of cases presented by Zhou et al. (2002) and the recovery imbibition curves of Morrow and Xie (2001) resulted from weak water-wet rocks.
In Figs. 2, 3, 4, 5 and 6, the data of the recovery factor are plotted in terms of logarithmic time, and the semi-logarithmic graphs are shown, respectively based on different datasets. Figure 2 shows recovery curves of weak water-wet rocks from the Morrow and Xie (2001) dataset. Figure 3 gives recovery curves resulting from the simulation of a single matrix block filled with oil and gas from the Ghasemi et al. (2020) dataset. Figure 4 shows recovery curves of strong water-wet rocks from the Standnes (2010a, b) dataset, Fig. 5 gives recovery curves of cases presented by the Fischer et al. (2006) dataset, and Fig. 6 shows recovery curves of cases presented by Zhou et al. (2002).

Results and discussions
To compare model accuracies, the average error values of all samples for each model are calculated based on the presented experimental data. Also, we divided sample timelines into three stages: the early stage, the middle stage, and the late stage. Then, we calculated the average error values of all samples for each stage. We called these errors local errors. This helps us to understand the model's advantages and disadvantages in predicting the recovery factor for each time stage. Also, RMSE was used to compare models with each other as it is a good parameter both for fitting purposes and model validation. The value of RMSE is more sensitive to local errors than the average errors that were defined.  Fig. 7, the value of local errors and average error for all models were compared. This graph shows us that the Logit-Hill and 1HMT models are more reliable ones in predicting recovery factors in a general way. Table 9 shows the values of RMSE and R 2 for the models used in different datasets. The best accuracy is attained when the value of parameter RMSE is close to 0, and the value of parameter R 2 is closer to 1. The highest value of R 2 was obtained by the Weibull model. Also, the lowest value of the RMSE parameter was calculated for this model (R 2 = 0.9992, RMSE = 0.0052). Table 10 shows the fitting coefficients of each model in different tests, which are shown from left to right, first twoparameter models and then three-parameter, respectively. "a" value can have negative values (a < 0) when one uses the two-parameter Probit model. Also, when both threeparameter models (1HMT and AHMT) are used, if the value of p in one model is equal to 1, it will have different values from 1 in the other model. Figures 8 and 9 show the values for RMSE and R 2 in models, respectively. The values of RMSE change in the approximate range of 0.0165-0.0270. In this range, the  lowest value was obtained for the 1HMT model, and the highest value was obtained for the Logit-Hill model. Also, for R 2 , the modeling values change in the approximate range of 0.988-0.995. In this range, the lowest value was obtained for the Logit-Hill, and the highest value was obtained for the 1HMT model. As it is known, the values of RMSE and R 2 for the 1HMT models are the minimum and maximum values, respectively. Thus, the best function for simulation among the five models is the 1HMT model. In Figs. 10 and 11, the results of the models fitted were compared to the experimental data for datasets 12 and 27. As can be seen from these graphs, and according to the calculations performed in the early, middle, and late times, the best models for fitting recovery factor curves in the spontaneous imbibition process are the 1HMT, Weibull, and Probit models, respectively.

Conclusions
In this research, the efficiency of different models for curve fitting of the recovery factor from fractured reservoirs during the spontaneous imbibition process was analyzed. The accuracy of these multi-parameter models was examined against various experimental data and using the coefficient of determination (R 2 ) and Root-Mean-Square Error (RMSE).
The most important results are as follows: • In contrast to previous studies that primarily tested the accuracy of one and two-parameter models, this research also explored the efficiency of three-parameter models.
The results showed that three-parameter models exhibit better accuracy in curve fitting of the imbibition data. The improved performance can be attributed to the presence of the p parameter (the third effective parameter). • The AHMT (All-Hit-Multi-Target), Logit-Hill, and Weibull models showed the best results for the early, middle, and late production times, respectively.
• The lowest average RMSE value was obtained for the 1HMT (One-Hit-Multi-Target) model, and therefore, this model is recommended for curve fitting the spontaneous imbibition processes in naturally fractured reservoirs.

Funding
The authors did not receive support from any organization for the submitted work.
Data availability All data generated or analyzed during this study are included in this published article. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.