Modeling of Brine/CO₂/Mineral Wettability Using Gene Expression Programming (GEP): Application to Carbon Geo-Sequestration

Abstract

Carbon geo-sequestration (CGS), as a well-known procedure, is employed to reduce/store greenhouse gases. Wettability behavior is one of the important parameters in the geological CO₂ sequestration process. Few models have been reported for characterizing the contact angle of the brine/CO₂/mineral system at different environmental conditions. In this study, a smart machine learning model, namely Gene Expression Programming (GEP), was implemented to model the wettability behavior in a ternary system of CO₂, brine, and mineral under different operating conditions, including salinity, pressure, and temperature. The presented models provided an accurate estimation for the receding, static, and advancing contact angles of brine/CO₂ on various minerals, such as calcite, feldspar, mica, and quartz. A total of 630 experimental data points were utilized for establishing the correlations. Both statistical evaluation and graphical analyses were performed to show the reliability and performance of the developed models. The results showed that the implemented GEP model accurately predicted the wettability behavior under various operating conditions and a few data points were detected as probably doubtful. The average absolute percent relative error (AAPRE) of the models proposed for calcite, feldspar, mica, and quartz were obtained as 5.66%, 1.56%, 14.44%, and 13.93%, respectively, which confirm the accurate performance of the GEP algorithm. Finally, the investigation of sensitivity analysis indicated that salinity and pressure had the utmost influence on contact angles of brine/CO₂ on a range of different minerals. In addition, the effect of the accurate estimation of wettability on CO₂ column height for CO₂ sequestration was illustrated. According to the impact of wettability on the residual and structural trapping mechanisms during the geo-sequestration of the carbon process, the outcomes of the GEP model can be beneficial for the precise prediction of the capacity of these mechanisms.

1. Introduction

Energy demands and industrial activities increase the amount of CO₂ emissions, the most specific greenhouse gas, into the atmosphere. To date, different methods, such as using biofuels as an alternative for fossil fuels, electrical energy as a clean source, and geological CO₂ sequestration, have been considered to limit the production and emission of CO₂. Carbon geo-sequestration (CGS), as a well-known procedure, is employed to reduce/store greenhouse gases [1]. Carbon geo-sequestration is contemplated as a relatively new method by which carbon dioxide could be mitigated effectively. Therefore, using this method’s untapped potential to reduce carbon emission to the atmosphere is highly recommended by many scientists [1,2]. Utilizing a comprehensive numerical method, Chen et al. [3] studied carbon capture and storage and its application in enhanced water recovery (CO₂-EWR) applications. Aiming to produce more brine while storing a higher amount of CO₂ underground, they applied different scenarios of CO₂ injection. They concluded that scenarios consisted in CO₂ enhanced water recovery in which the co-injection of brine and pre-injection of brine were used had a better performance due to the fact that not only had more CO₂ been stored, but pressure changes could also be controlled more effectively. Furthermore, it was found that when more injection wells were used, a larger amount of CO₂ could be stored underground. However, it was concluded that drilling more injection wells will significantly increase the costs of carbon capture and storage (CCS) practice, and an implementation of this strategy is not cost-effective. The possibility of the co-injection of impurities along with CO₂ allows for the direct disposal of flue gas, and hence a significant reduction in the cost of CO₂ sequestration projects by eliminating the separation process. Based on this, different studies have sought to examine the feasibility of the sequestration of CO₂-N₂ or CO₂-SO₂ mixtures in saline aquifers [4,5].

The effect of dip angle and salinity of brine on the amount of stored CO₂ has been investigated by Jing et al. [6]. In order to simulate the process, a three-dimensional model was developed. It was found that as the salinity of underground water increased, the amount of stored CO₂ decreased. Regarding the dip angle, it was found that as the dip angle increased, CO₂ migration distance surged. To maximize the performance of a CO₂ geo-sequestration practice, a reservoir with a smaller dip angle and lower salinity should be given priority in the selection process.

During CO₂ geo-sequestration, natural and artificially introduced fractures could impose a risk to the success of the operation on the grounds that CO₂ could migrate upward and return to the atmosphere. In order to control CO₂ migration and trap the injected CO₂ forever, a wide range of physicochemical mechanisms have been investigated, including structural trapping [7], mineral trapping [8], dissolution trapping [9], and residual trapping [10]. Physical adsorption methods could also be used in CO₂ storage in sandstones. Using this approach, carbon dioxide can be stored in the interior layers or on the surface of clays present in sandstone formations [11,12]. To ensure secure carbon sequestration, the selection of appropriate geological formations must be made carefully. Otherwise, carbon dioxide may not be trapped effectively.

Regarding the performance of various mechanisms, it has been found that a profusion of different factors, such as aquifer and cap rock properties, are of the utmost importance to be considered carefully. Additionally, a widely held belief is that trapping mechanisms could be either active or inactive in different periods. At the early stages of a sequestration operation, not only are structural and residual trapping mechanisms salient, but it is also believed that they are the only affecting mechanisms [13,14]. These mechanisms are illustrated in Figure 1. In structural trapping, carbon dioxide is trapped below a seal layer with an extremely low permeability, and residual trapping is relevant to the cases in which there is no cap rock [15]. Other mechanisms for the storage capacity of CO₂ sequestered in oil reservoirs is mineral trapping and solubility. The solubility of CO₂ in remaining oil is much higher than that of formation water [16]. It should be noted that the heterogeneity of rock permeability can strongly affect all trapping mechanisms, the details of which can be found elsewhere [17]. In the residual trapping mechanism, CO₂ is not mobile due to capillary forces, which are highly dependent on certain salient influencing factors, such as the initial saturation of CO₂, the morphology of reservoir rock, the interfacial tension between CO₂ and brine, and the wettability of the reservoir rock. As a result, the role of wettability in CO₂ sequestration capacity and prediction of CO₂ leakage to the atmosphere has been investigated in the past [15,18,19]. The crucial role of wettability in the mobility of different phases is one of the important research topics in petroleum engineering, and its alteration according to any introduced variation in reservoir conditions should be considered carefully.

Several studies have dealt with wettability behavior by determining the contact angle between brine and carbon dioxide on various minerals under different conditions, such as pressure, the salinity of the brine, and temperature. Farokhpoor et al. [21] determined the contact angle of a CO₂/brine system on a range of different minerals, such as calcite, feldspar, mica, and quartz. They found that, while the variation of pressure influenced the wettability of mica and changed it from a strong water-wet material to a weaker one, it had no impact on other minerals. In addition to that, it was found that some minerals, such as feldspar, quartz, and calcite, had a maximum contact angle of 36 degrees at critical pressure. Chen et al. [22] studied the contact angle of water on silica at various temperatures ranging from 318 to 383 K and pressures between 2.8 and 32.6 MPa using the molecular dynamic simulation method. It was found that to control and modify the effect of pressure and temperature on the wettability of minerals, surface functional groups could be employed successfully. They investigated various contact angles of the CO₂/brine system on quartz at different pressures, salinity ranges, and temperatures. Wettability was simulated by the implementation of a molecular dynamic simulation. It was concluded that while contact angle is not a strong function of temperature and pressure, ionic strength has a direct influence on the contact angle of water. A comparison between simulation and obtained experimental contact angles showed that the simulation predicted the contact angles precisely. In another work, in order to calculate the interfacial tension in a CO₂-brine system, an empirical method was employed by Mutailipu [23]. Contact angles of supercritical, liquid, and gaseous CO₂ on different minerals, such as limestone, quartz, and Brea sandstone, were calculated in terms of salinity, temperature, and pressure. It was found that limestone and Berea sandstone experienced alteration to less water wet rocks, whereas quartz remained relatively unchanged when supercritical conditions were dominant.

Based on the above studies, investigating the wettability of carbon dioxide can be considered a contributing and paramount pathway for other related studies and investigations. Hence, the amount of CO₂ which could be stored in a reservoir by structural trapping and the residual mechanisms could be predicted. Unfortunately, it was found that there are a few comprehensive models in the literature by which the brine/CO₂/mineral contact angle could be calculated at different reservoir conditions. Additionally, large differences and uncertainties in the studies have made it challenging to develop an appropriate model to predict wettability in the above-mentioned system. This involves the development of a precise and comprehensive model by which the CO₂ contact angle can be calculated. In recent decades, new models have been developed to solve complicated systems. Using soft computing methods, the profusion of different problems has been solved and many outstanding methods have been proposed [24,25,26,27,28,29]. Although these smart models are very useful in solving complicated systems, they typically encounter some limitations which are intrinsic to them, such as over/under fitting problems and the existence of a “black box”, which is a necessity of these methods.

In order to predict the contact angle in ternary systems of brine, CO₂, and minerals, intelligent models have been proposed, which are black box and need specific software, such as Matlab or python for calculations [20,30]. The present study deals with some of the above-mentioned limitations using a highly valuable method known as gene expression programming (GEP) to estimate the wettability in a brine/CO₂ system on various minerals, such as feldspar, mica, quartz, and calcite. It is worth mentioning that the use of GEP method in this study in order to develop accurate models for predicting the wettability of brine/CO₂/mineral is due to its advantages, mainly in terms of its accuracy and ability to generate explicit and user-friendly correlations that can be integrated in other applications. The wettability of various minerals regarding the brine and CO₂ system has a great importance in carbon dioxide geological storage since it can considerably affect the residual and structural trapping. Consequently, simulating the wettability behavior of the brine/CO₂/mineral system is vital when operating conditions have an influence on this system. Figure 2 is a general sketch of the methodology used in the present research. To this end, a large dataset comprising 630 values of contact angles in various reservoir conditions was gathered from the literature. Furthermore, the untapped potential of leverage methods has been utilized to ensure the validity of the proposed model by considering the influence of different variables on wettability, including pressure, temperature, and salinity. The main contribution and novelty of this study consists of establishing a user-friendly correlation for predicting contact angle of brine/CO₂/mineral in a ternary system under extensive operational conditions. To the best of our knowledge, no previous work has implemented the GEP technique for predicting the wettability behavior of the ternary system of brine/CO₂/mineral.

The remaining sections of this study are outlined as follows: First, the theoretical background involving data collection and the principles of the applied soft computing approach, namely gene expression programming, is presented. Then, the implementation procedure is highlighted, and the results of the effect of operational parameters on contact angles, applicability domain, and sensitivity analysis are presented and discussed. Next, the application of the proposed model for calculating the CO₂ column height in the subsurface is illustrated. The study ends with conclusions, which recap the main results of the investigation.

2. Theoretical Background

2.1. Data Collection

To develop a representative model for the brine/CO₂/minerals system, we utilized an extensive dataset comprising 630 values of contact angles at various reservoir conditions. Recently, numerous attempts have been made to investigate this issue and many contact angles in receding, static, and advancing stages have been reported. Many of the contact angles are unreliable because the values obtained were significantly different although the conditions were similar. The differences in measurements reported could be a result of surface contamination by other minerals or even the roughness of the surface [31,32]. Contaminants on the surface could be a major source of the problem. As a result, a wide range of genuine approaches have been invented and utilized to clean the surface of minerals from contamination. Problems which are due to surface roughness should also be addressed properly. Hysteresis, the surge of contact angle on hydrophobic surfaces and its dip on hydrophilic surfaces, also represents a major problem caused by increased surface roughness [33,34].

Using the above-mentioned studies [31,32,33,34], a large set of data points was gathered from the literature, which is utilized in the current study. In order to make calculations more precise and take previously mentioned factors into account, we used the factor (θ⁰) proposed by Daryasafar et al. [20]. Theta zero represents the wettability of minerals when considering the impact of surface contamination and roughness. This factor is expressed as below [20]:

$${\mathsf{\theta }}^0=\mathrm{round}\left(\frac{{\theta }_i }{10}\right)$$

(1)

where ${\theta }_i$ indicates the contact angle between brine, mineral, and carbon dioxide in an environment with a salinity equal to zero at ambient temperature and pressure. As expected, a direct relationship between the proposed factor and the system’s contact angles is evident. Additionally, knowing that the contact angle is highly dependent on surface features, we believe that theta zero, which shows the impact of surface roughness and contaminants, must be studied carefully.

The present study sheds light on the development of accurate paradigms for estimating contact angles. In order to develop the models, valuable information about salinity (M), minerals (calcite, quartz, mica, and feldspar), temperature (K), pressure (MPa), type of contact angles, and the defined factor (θ⁰) was collected and introduced to the presented model as input variables, similar to the work of Daryasafar et al. [20]. As described below, the contact angle is a function of the above-mentioned variables:

$$\mathsf{\theta }\left(\mathrm{adv} \mathrm{or} \mathrm{rec} \mathrm{or} \mathrm{st}\right)=\mathrm{f} \left(\mathrm{P},\mathrm{T},\mathrm{mineral} \mathrm{type}, \mathrm{salinity},{ \mathsf{\theta }}^0, \mathrm{contact} \mathrm{angle} \mathrm{type}\right)$$

(2)

Various types of contact angle can be introduced to the developed model using 1 for static contact angle, 2 for advancing, and 3 for receding type. Additional information representing the data bank is given in

Table 1. Statistical description of the employed database.

	Pressure (MPa)	Temperature (K)	Salinity (M)	Number (sta/adv/rec)	Theta. zero	Contact Angle
Minimum	0.04	296	0	1	0	6.84
Average	9.59	320.52	1.27	1.95	1.73	36.85
Maximum	40.05	373	12.88	3	4	122.32
Median	9.02	318	0.2	2	2	32.59
Mode	10	308	0	1	3	68.03
Kurtosis	3.05	1.09	7.73	−1.74	−1.22	0.96
Skewness	1.32	1.17	2.65	0.08	−0.16	0.96

.

It is worth mentioning that to ensure the highest degree of accuracy and robustness in the generated GEP-based correlations, the amassed database was split into four groups according to mineralogy, namely calcite, feldspar, mica, and quartz. Besides, it is necessary to add that for each mineralogy group, another splitting of the data was performed randomly into training (80% of the data) and testing (20% of the data) sets.

Figure 3 indicates the box plots of the input parameters, including pressure, temperature, salinity, contact angle type and theta zero. The number of input data for each mineral involve 73, 30, 145, and 382 for calcite, feldspar, mica, and quartz, respectively. These plots were drawn vertically to obtain a better comparison between the data set groups. For the pressure input (Figure 3a), the median values of all minerals are close to each other, and feldspar has a wide variation domain without any outlier data, while the others possess one or more outlier data. For the temperature parameter (Figure 3b), the box height of quartz mineral data indicates the largest distribution in comparison to the other minerals, especially feldspar. Both quartz and feldspar have one outlier datum. The most outlier data was obtained for salinity input with the narrow size of the box heights for all minerals (Figure 3c). The box plots of the minerals for contact angle type input showed different distributions without any outlier data (Figure 3d). The results also showed that significant differences were found between feldspar and mica data for the theta zero parameter (Figure 3e). The middle quartile domain or box height for mica was large, which confirms that it is the most divergent toward other minerals.

2.2. Gene Expression Programming (GEP)

Gene expression programming (GEP), which is very useful in developing white-box models, could be used in expanding computer programs. This evolutionary-based and advanced approach is an appropriate method by which different systems could be described using inputs and desirable outputs. The GEP was first developed and introduced by Ferrera [35,36] and has been contemplated as a new variation of genetic programming [37]. In the older version of GEP, problems, such as faulty explorations and a limited number of regression methods, have been addressed [35,36]. Regarding the fundamentals of this method, the expression tree (ET) and the chromosome should be deemed as the GEP’s conceptualization and data processing basis [38,39]. Utilizing previously obtained data from experimental works or simulations, the GEP method can derive real solutions (i.e., ET) from the chromosomes. As illustrated in Figure 4, the genes comprise a fixed-length symbolic inventory in which mathematical operators and terminal variables are depicted [40].

In almost all genetic strategies, it is of the utmost importance to initialize the program using a logical process by which the GEP will start moving towards a solution. Employing a fitness function and genetic operators, the GEP can evaluate the newly created genes and, resultantly, it could find appropriate candidates. To complete the task, the GEP implements a wide range of genetic operators, e.g., tournament-based selections, mutation operators, and crossover. Furthermore, transposition and recombination operators are also worthwhile, and all operators will be reused endlessly and iteratively unless the cease criterion is met. Figure 5 represents the main steps for designing a GEP structure to obtain the best model.

2.3. Model Development

In the present work, the selected experimental database was randomly divided into two parts, namely training and testing data sets, for developing the GEP model. Hence, 80% of the total data was used as the train set, while the remaining 20% was considered as test set. In addition,

Table 2. The utilized setting parameters for the implemented correlations in this study.

Parameters	Value/Setting
Population size	60
Crossover’s probability	90%
Mutation’s probability	15%
Elitism	10%
Type of selection	Linear ranking
Max. number of generations	100

represents the considered key parameters of the applied evolutionary algorithm.

To evaluate the accurateness and performance of the improved models, certain statistical parameters were utilized consisting of average percent relative error (APRE), average absolute percent relative error (AAPRE), root mean square error (RMSE), standard deviation error (STD), and coefficient of determination (R²). Definitions and equations of those parameters are given below:

A.

Average percent relative error (APRE). It measures the relative deviation from the experimental data and is defined by:

$$\%APRE=\frac{100}{M}{\displaystyle \sum }_{i=1}^M\frac{\left({\theta }_{i_{exp}}-{\theta }_{i_{pre}}\right)}{{\theta }_{i_{exp}}}$$

(3)

B.

Average absolute percent relative error (AAPRE). It measures the relative absolute deviation from the experimental data and is defined as:

$$\%AAPRE=\frac{100}{M}{\displaystyle \sum }_{i=1}^M\left|\frac{\left({\theta }_{i_{exp}}-{\theta }_{i_{pre}}\right)}{{\theta }_{i_{exp}}}\right|$$

(4)

C.

Root mean square error (RMSE). It measures the data dispersion around the zero deviation and is defined by:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^M{\left({\theta }_{i_{exp}}-{\theta }_{i_{pre}}\right)}^2}{M}}$$

(5)

D.

Standard deviation (SD). It is a measure of dispersion, and a lower value shows a smaller degree of scattering. It is defined as:

$$STD=\sqrt{\frac{1}{M-1}{\displaystyle \sum }_{i=1}^M{(\frac{\left({\theta }_{i_{exp}}-{\theta }_{i_{pre}}\right)}{{\theta }_{i_{exp}}})}^2}$$

(6)

E.

Coefficient of determination (R²). It is a simple statistical parameter that exhibits how a good model matches the data. The closer the R² value is to 1 confirms the better fitting of the model. It is defined as:

$$R^2=1-\frac{{{\displaystyle \sum }}_i^M{({\theta }_{i_{exp}}-{\theta }_{i_{pre}})}^2}{{{\displaystyle \sum }}_i^M{({\theta }_{i_{pre}}-\overline{{\theta }_{i_{exp}}})}^2}$$

(7)

In the above equations, $\theta$ points out the wettability values, the subscripts exp and pre denote the experimental and predicted values, $\overline{\theta }$ is the average value of the wettability, and M is the number of data points.

3. Results and Discussion

3.1. Model Implementation

In this study, the performance of the implemented correlations, the GEP model, was evaluated to predict the contact angle (advancing, receding, and static) of brine/carbon dioxide/rocks under various operating conditions, including pressure, temperature, and salinity. To establish an accurate correlation to predict the contact angle in the brine/CO₂/minerals ternary system using the GEP approach, the control parameters were appropriately tuned. In addition, the generated correlations for all four minerals are listed in

Table 3. Final explicit correlations for estimating the brine/CO₂/mineral wettability.

Mineral/Parameters	Correlations
Calcite
$\theta$	$A_1\times P+A_2+A_3+A_4+472.2$
$A_1$	$56.28-0.09045\times \left(T+{\theta }_o\right)-0.006978\times P^2\times \mathit{ln}\left({\theta }_o\right)-4.365\times N\times \mathit{tanh}\left({\theta }_o\right)$
$A_2$	$\left(\frac{N}{P}\right)\times \left(\frac{9.537\times {10}^{-5}\times N^2}{P^2}-5.1\right)-\frac{25.01\times P}{N}+\frac{2314}{T}+\frac{0.01692\times T\times N}{\mathit{tanh}\left(P\right)}$
$A_3$	$-2.808\times S-17.13\times N-2.292\times {\theta }_o\times \left(1-S\right)$
$A_4$	$-24.68\times \mathit{tanh}\left(S\times N\right)-78.51\times \mathit{ln}\left(T\right)+22.29\times \left(\mathit{tanh}\left(S\right)+\mathit{ln}\left({\theta }_o\right)\right)$
Feldspar
$\theta$	$\frac{A_1}{P}+A_2+A_3+A_4-16.72$
$A_1$	$2.434\times \sqrt{T}+19.51\times \mathit{ln}\left(P\right)-37.74\times \mathit{tanh}\left(P\right)$
$A_2$	$\frac{6.679\times \sqrt{T}}{4.504\times P-35.67}-\frac{7.525\times P^2}{18.01\times P-142.4}$
$A_3$	$\frac{462.3\times P}{T}-\frac{1.6683\times P^2}{P+S}$
$A_4$	$0.05923\times T\times \mathit{tanh}\left(P\right)+4.81\times \sqrt{P}+0.0004529\times P\times T\times \mathit{exp}\left(S\right)$
Mica
$\theta$	$A_1\times {\theta }_0^3\times \left(N-{\theta }_0\right)+A_2\times P+A_3+A_4+A_5+2554$
$A_1$	$1.482-0.04563\times P$
$A_2$	$664.6-68.71\times \sqrt{T}-0.02321\times P\times {\theta }_0^2$
$A_3$	$1.807\times \left(S+T\times P\right)+0.9754\times N\times \left(N+S\times {\theta }_0\right)+1.234\times N\times {\theta }_0\times \left(N+P\right)$
$A_4$	$34.58\times \sqrt{P}-12.37\times \mathit{tanh}\left(P^2\times S\right)-2564\times \mathit{tanh}\left(P+N\right)$
$A_5$	$-\frac{15.56\times {\theta }_0\times \sqrt{P}\times \mathit{ln}\left(P\right)}{N}-\frac{4959\times \left(P+N\right)}{T}$
Quartz (T > 300 K)
$\theta$	$A_1\times \sqrt{T}+A_2\times \sqrt{S}+A_3\times T+A_4\times P+A_5+A_6+3749$
$A_1$	$53.05\times P-3.406\times \sqrt{S}\times \sqrt{{\theta }_0}-611.1\times \sqrt{N}$
$A_2$	$40.25\times {\theta }_0+12.93\times N-23.62$
$A_3$	$34.26\times \sqrt{N}-5.966\times N-11.69$
$A_4$	$1909-496.1\times \mathit{ln}\left(T\right)-0.02984\times P\times S\times N$
$A_5$	$0.2824\times S\times \left(P-{\theta }_0\right)\times \mathit{ln}\left(P\right)-0.5521\times {\theta }_0\times \left(P+{\theta }_0\right)\times \mathit{tanh}\left(N\right)$
$A_6$	$1886\times N+\frac{1712.597\times \left(P+{\theta }_0\right)}{\left(P+T\right)}$
Quartz (T < 300 K)
$\theta$	$A_1\times {\theta }_0\times N^2+A_2+A_3+A_4+276.2$
$A_1$	$P\times N\times \left(0.09477\times N-0.2665\times \mathit{tanh}\left(P\right)\right)-0.7928\times {\theta }_0\times \mathit{tanh}\left(N\right)$
$A_2$	$8.188\times P-0.7642\times \left(T+\mathit{ln}\left(P\right)\right)+3.523\times \left(S+N\times \left(1-S\right)\right)$
$A_3$	$18.22\times ({\theta }_0-\mathit{ln}\left(T\right))-\frac{7.106\times P}{N}-1.082\times P\times N^{\frac{3}{2}}$
$A_4$	$49.59\times \mathit{tanh}\left(N\right)+15.23\times \mathit{tanh}\left(6.112\times T\times S\right)-11.87\times \mathit{tanh}\left(P\times S\times N\right)+15.03\times \mathit{tanh}\left(P\times S\times {\theta }_0\right)$

.

The contact angle prediction capability of the proposed model was confirmed by preparing different measurements indexes, such as APRE, AAPRE, RMSE, STD, and R². Both statistical evaluation and graphical analyses were performed to show the correctness and performance of the developed model. The statistical parameters are listed in

Table 4. Statistical evaluation of the implemented correlations.

		Calcite	Feldspar	Mica	Quartz (T > 300 K)	Quartz (T < 300 K)
Training set	APRE (%)	−0.5	0.71	14.71	23.62	13.45
	AAPRE (%)	4.4	1.69	14.71	27.54	14.28
	RMSE	1.81	0.37	7.340	10.97	3.2
	STD	0.004	0.0004	0.026	0.110	0.037
	R2	0.996	0.972	0.983	0.806	0.951
Test set	APRE (%)	−7.12	−0.001	13.38	21.11	12.54
	AAPRE (%)	10.55	1.01	13.38	27.71	12.54
	RMSE	3.29	0.23	7.020	10.49	3.180
	STD	0.04	0.0002	0.022	0.112	0.020
	R2	0.985	0.974	0.965	0.678	0.978
Total	APRE (%)	−1.86	0.57	14.44	23.12	13.26
	AAPRE (%)	5.66	1.56	14.44	27.57	13.93
	RMSE	2.2	0.340	7.280	10.87	3.190
	STD	0.013	0.0004	0.025	0.110	0.034
	R2	0.994	0.972	0.980	0.780	0.956

for calcite, feldspar, mica, and quartz (below 300 K and over 300 K). A graphical comparison between different minerals in predicting contact angle is shown in Figure 6. Table 4 and Figure 6 demonstrate that the results estimated by the presented models are in good accordance with the experimental data. It is evident from the statistical evaluation that the GEP correlations have excellent prediction performance with overall AARE% values of 5.66%, 1.56%, 14.44%, and 13.93%, for calcite, feldspar, mica, and quartz, respectively.

A clear explanation for the accuracy and proficiency of intelligent models can be achieved through visual comparison between predicted and actual experimental data. For this purpose, the anticipated contact angles of the brine/carbon dioxide system on various mineral surfaces were plotted versus actual values for both train and test subsets, as shown in Figure 7. This depiction indicates that the cloud of estimated values of brine/CO₂ contact angle by the GEP model is accumulated around the bisector line of y = x, which confirms the acceptable accuracy of the proposed model. Based on this fact, the developed GEP model is accurate and has an excellent capability for simulating the behavior of brine/CO₂/mineral wettability. For additional confirmation, the relative deviation of predicted data was plotted versus the experimental contact angles data. As shown in Figure 8, the low ranges of distribution around the zero-line error express the good capacity and proficiency of GEP as an estimative tool.

According to the above-mentioned statistical analyses, the GEP models can show significant compatibility between the estimated and experimental values since there is excellent overlap between the predicted and target data points. To establish a better description of the absolute relative errors of the implemented model, the depiction of cumulative data frequency against absolute percent relative error is drawn in Figure 9. This plot is described as the percent of total data points used and illustrates that more than 80% of the estimated contact angles have an absolute relative error less than 25%.

For further verification of the proposed models’ validity and investigation of the average absolute relative error over various input parameters, Figure 10 is presented. According to this figure, the proposed GEP model can accurately estimate the experimental data within different ranges of input parameters. This indicates the robustness and smartness of the presented model for broad domains of operating conditions. The amount of AAPRE is less than 22% for almost all ranges of input variables, except pressure over 30 MPs and theta zero = 1, thus confirming the above statement. On the other hand, choosing a reliable model is one of the main issues for estimating any parameter within a specific range of input variables. Figure 10 demonstrates the capability of the GEP model in predicting the wettability of brine/CO₂/rocks for various ranges of operating conditions with high precision.

3.2. Effect of Operational Parameters on Contact Angles

In this study, many attempts have been made out to show the ability and accuracy of the proposed GEP model in trend to estimate brine/CO₂ contact angles for different minerals under wide rages of pressure, temperature, and salinity. Figure 11 illustrates a comparison of the contact angle with experimental data by considering the variation of influencing parameters in the proposed models. As can be observed in Figure 11a, the measured values for contact angles versus pressure indicate that the implemented model has high accuracy in the prediction of the process trend. In addition, the depiction demonstrates that the brine/CO₂ contact angle decreases with increasing pressure. The influence of temperature on wettability behavior was investigated as well and the estimation ability of the developed model was confirmed based on Figure 11b. The GEP model shows sufficient precision in estimating brine/CO₂ contact angle, where it decreases with temperature up to 340 K and increases at high temperatures.

Salinity is another operating parameter that affects the performance of the proposed model in the trend estimation of brine/CO₂ contact angle. According to the used datasets, the GEP model is excellent for the trend predictions of wettability behavior versus salinity. It is evident from Figure 11c that contact angles increase by enhancing NaCl concentration. Results obtained by other researchers confirm this trend [20]. However, different results for these trends are reported by researchers, e.g., for the effect of temperature on brine/CO₂/quartz [41]. Based on the results of the developed intelligent model, one can conclude that the GEP has good performance and shows reliable and accurate results.

3.3. Applicability Domain and Sensitivity Analysis

Some outliers are typically included as an intrinsic feature of almost any dataset. Outliers show a behavior different from the rest of the data points in a dataset. Believing that the model will return more precise and reliable results provided that the outliers are detected, we have implemented the leverage algorithm in the current study. In this approach, the deviation of model outputs from its corresponding experimental data is named as standardized residuals. In addition, Hat matrix is described as follows [42,43,44]:

$$h=X{\left(X^{T}X\right)}^{-1}{X}^T$$

(8)

In this equation, X represents the $p\times q$ matrix (in which p and q stand for the number of actual data points and dimension of the model, respectively) and $X^T$ represents the transpose matrix. Hat vector is determined by the diagonal elements of the Hat matrix. In addition, the leverage limit (h^∗), as a warning value, is defined by [45]:

$$h^{*}=\frac{3\left(q+1\right)}{p}$$

(9)

Clarification of the doubtful data points was investigated by a Williams plot, in which standardized residuals are plotted versus hat values, as shown in Figure 12. The leverage limits (h^*) corresponded to 0.25, 0.4, 0.125, and 0.047 for calcite, feldspar, mica, and quartz, respectively. Six data points had hat values above these thresholds. It is evident from Figure 12 that just five outlier data points were laid out of the reliable zone (outside the range −3 to 3) for all minerals. Although these data points appear differently from others, they did not show a strong effect on the performance of the proposed model. Therefore, all the experimental data were used.

To study the sensitivity analysis, the so-called Pearson correlation coefficient $r_i$ was utilized to show the influence of input parameters on the prediction of the contact angle of the brine/CO₂/minerals system. This relevancy factor ($r_i$) is defined for each input parameter, i, as follows [46,47]:

$$r_i=\frac{{{\displaystyle \sum }}_{k=1}^M\left(A_{i,k}-{\overline{A}}_i\right)(B_k-\overline{B})}{\sqrt{{{\displaystyle \sum }}_{k=1}^M{\left(A_{i,k}-{\overline{A}}_i\right)}^2{{\displaystyle \sum }}_{k=1}^M{\left(B_k-\overline{B}\right)}^2}} , \left(i=1,\dots ,5\right)$$

(10)

where $A_{i,k}$,$\overline{A}$,$N$, $B_k$, and $\overline{B}$ represent input parameters, average of inputs, number of the data points, output parameter, and average of output, respectively. The value of $r_i$ is located within −1 to 1 and the large values correspond to the strong correlation. Moreover, the increasing or decreasing output parameter with variations in $M_i$ attribute to a positive or negative sign, respectively.

Figure 13 shows the identified relevancy factors for all five input parameters, which affect the amount of wettability in CO₂, saline water, and minerals. All these parameters have a significant contribution to the value of the contact angle. Temperature possessed the least amount of influence compared to the others, negatively (r~0.035). The r-value of salinity shows the utmost significant dependence (r = 0.23), followed by pressure (r = 0.183). Due to the obtained values for relevancy values, it can be concluded that salinity and pressure are the utmost impressive parameters, and temperature is the least significant for identifying the amount of contact angle in the brine/CO₂/mineral system.

3.4. Application of the Proposed Model for CO₂ Sequestration

In the context of showing the impact of estimating wettability using the proposed GEP-based correlation on the trapping mechanisms, we performed a comparison between the different storage heights associated with the structural trapping capacity in mica mineral type. The illustration of the trapping mechanism can be found elsewhere [27,48]. The following expression was used for calculating the CO₂ column heights in underground conditions [27,49]:

$$H=\frac{2\gamma \cos \left(\theta \right)}{∆\rho gR}$$

(11)

where g is the gravitational acceleration, R denotes the pore throat radius, $\gamma$ is the interfacial tension (IFT) of the system CO₂–brine, $\theta$ points out the contact angle between CO₂, brine, and the rock surface, and $∆\rho$ represents the density difference between CO₂ and brine.

The compared storage heights included the real values of height calculated by considering the wettability measurements reported in the study of Arif et al. [49], the values of height calculated by considering the estimations of wettability using the suggested GEP-based correlation, and the values of height calculated by neglecting wettability (i.e., $\theta =0$°). The calculations were carried out at two different temperatures (323 K and 343 K) and various pressure values. The IFT values considered in Equation (11) under the associated conditions of pressure and temperature were obtained from Arif et al. [49]. It is worth mentioning that for the calculation of wettability using GEP-based correlations for the considered conditions, we assumed the following parameters: ${\theta }_o=0$, receding contact angle ($N$ = 3), and Salinity = 3.42 M (20 wt% NaCl brine [49]). The results of the comparison are shown in

Table 5. Comparison between CO₂ column height values of the considered cases.

Consideration of the wettability measurements reported in the study of Arif et al. [46]	Temperature (K)	Pressure (MPa)	$∆\rho$(kg/m3)	${\displaystyle \mathrm{IFT}, \gamma }$(mN/m)	$\theta$(°)	CO2 column height (m)
	323	5	1031	55	29	952
		10	755	43	50	747
		15	445	38	67	681
		20	359	36	74	562
		25	320	33	79	402
	343	5	1032	58	26	1031
		10	881	45	41	786
		15	625	40	58	691
		20	474	38	68	613
		25	380	36	74	533
	Temperature (K)	Pressure (MPa)	$∆\rho$(kg/m3)	${\displaystyle \mathrm{IFT}, \gamma }$(mN/m)	$\mathrm{GEP} \mathrm{estimated} \theta$(°)	CO2 column height (m)
Consideration of the estimations of wettability using the suggested GEP-based correlation	323	5	1031	55	14.04	1056.18
		10	755	43	36	940.34
		15	445	38	51	1096.73
		20	359	36	61	992.16
		25	320	33	70	719.81
	343	5	1032	58	13.62	1114.72
		10	881	45	32	884.02
		15	625	40	44	939.55
		20	474	38	52	1007.28
		25	380	36	57	1053.01
	Temperature (K)	Pressure (MPa)	$∆\rho$(kg/m3)	${\displaystyle \mathrm{IFT}, \gamma }$(mN/m)	$\theta$(°)	CO2 column height (m)
$\mathrm{Neglecting} \mathrm{wettability} (\mathrm{i}.\mathrm{e}., \theta =0$°)	323	5	1031	55	0	1088.699
		10	755	43	0	1162.319
		15	445	38	0	1742.72
		20	359	36	0	2046.501
		25	320	33	0	2104.592
	343	5	1032	58	0	1146.97
		10	881	45	0	1042.41
		15	625	40	0	1306.12
		20	474	38	0	1636.10
		25	380	36	0	1933.40

and Figure 14. As can be seen, the estimation of brine/CO₂/mineral wettability using the GEP-based correlations yields acceptable values for CO₂ column height compared with the case when brine/CO₂/mineral wettability was neglected. It is necessary to add that in many published studies [27], this vital parameter was neglected in the calculation of CO₂ column height, and this has a significant effect on this calculation as reported in Table 5 and Figure 14. Therefore, the newly proposed correlations for predicting brine/CO₂/mineral wettability in this study are of great interest to several calculations related to carbon geo-sequestration.

Lastly, it should be highlighted that the newly proposed GEP-based correlations for the modeling of brine/CO₂/mineral wettability are recommended for cases that fall within the ranges of application. Indeed, these explicit correlations can be applied for cases described by conditions that are outside of this applicability realm, but with careful attention as its accuracy can vary from one case to another. Nevertheless, as an extensive database with wide-ranged conditions was involved in the development of the correlations, these correlations can be invoked for predicting the wettability of many brine/CO₂/mineral systems having proprieties located within the range of the input parameters mentioned above.

4. Summary and Conclusions

In the present study, a machine learning algorithm, the so-called GEP technique, was developed to predict contact angles of carbon dioxide and saline water on different minerals under a wide range of environmental conditions. Pressure, temperature, salinity, mineral type (calcite, quartz, mica, and feldspar), type of contact angles, and the defined theta zero factor were used as input parameters and the contact angles of brine/CO₂ constituted the model output. The obtained values from the GEP model were compared with actual experimental data and the following important conclusions can be drawn based on the achieved outputs:

• Different measurements indexes, such as APRE, AAPRE, RMSE, STD, and R², confirmed the reliability and accuracy of the implemented model.
• Average absolute percent relative errors of the implemented model proposed for calcite, feldspar, mica, and quartz were obtained 5.66%, 1.56%, 14.44%, and 13.93%, respectively, which confirms the significant performance of the GEP algorithm.
• The GEP correlation was able to predict more than 80% of the considered data points with ARE less than 25%.
• The applied data points did not show significant outliers, and the proposed GEP model was successful in the trend estimation of brine/CO₂ contact angles for different minerals under wide ranges of pressure, temperature, and salinity.
• Investigation of sensitivity analysis indicated that the contact angles of brine/CO₂ on various minerals could be positively affected by salinity and pressure and negatively by temperature.
• According to the impact of wettability on the residual and structural trapping mechanisms during the carbon geo-sequestration process, the outcomes of the GEP model in this study can be beneficial for the precise prediction of these mechanisms’ capacity.