Enhanced water saturation estimation in hydrocarbon reservoirs using machine learning

Data collection & processing

In this study, the data required for developing high-accuracy and efficient ML models were extracted from hydrocarbon fields located in the southwestern region of Iran. The data acquisition process focused on obtaining high-quality, real-field measurements to ensure the reliability and representativeness of the dataset. The dataset includes a broad spectrum of well-log parameters such as SW, DEPT, HTNP, RT, CGR, SGR, HCAL, DTCO, RHOZ, and TEMP, comprising a total of 30,660 individual data points. This comprehensive dataset enables the analysis of variable behavior under different reservoir conditions and allows ML models to effectively predict and evaluate complex petrophysical relationships. Table 1 provides a detailed description of the input features used in this study.

The dataset used in this study was compiled from various published sources and internal databases to ensure diversity and generalizability of the machine learning models. Although the primary focus was on petrophysical log data, the geological setting of the wells reflects typical characteristics of sandstone-dominated hydrocarbon reservoirs. The reservoir formations include siliciclastic sequences with varying degrees of sorting and compaction, resulting in a wide porosity range from approximately 10 to 25%. These variations are influenced by depth, depositional environment, and cementation effects. The lithology primarily consists of medium to fine-grained sandstones interbedded with shales, which introduces heterogeneity in fluid saturation and log responses. Such geological complexity makes conventional water saturation estimation less reliable and further justifies the application of machine learning techniques to improve prediction accuracy in these heterogeneous formations.

Prior to model training, input features were standardized using z-score normalization to ensure consistent scaling across all variables. This process involved subtracting the mean and dividing by the standard deviation of each feature, resulting in standardized values with zero mean and unit variance. Standardization was especially important for models such as SVM and LR, which are sensitive to feature scale. Furthermore, Pearson correlation analysis was performed to identify potential multicollinearity among input features. Despite moderate correlations between some parameters, all nine features were retained due to their individual importance and non-redundant contributions to the target variable (Sw) prediction.

For data analysis, various tools such as scatter plots, histograms, heatmaps, box plots, and violin plots were utilized. Scatter plots are used to depict the relationship between two numerical variables and can assist in identifying correlations, clusters, and outliers. In contrast, histograms represent the distribution of a numerical variable by dividing the data into intervals and displaying the frequency of each interval as a bar chart. Heatmaps use colors to indicate the intensity or magnitude of data in a two-dimensional matrix, making them highly useful for analyzing correlation matrices and multidimensional data. Box plots summarize data distributions by presenting the median, quartiles, and interquartile range, which are helpful in identifying outliers and comparing distributions. Additionally, violin plots, which combine box plots and density curves, not only show the median and quartiles but also represent the shape of the data distribution, making them particularly useful for detailed analysis and identifying multimodal distributions.

In this paper, the outlier removal method was applied using a Gaussian-based approach to enhance data quality. To visualize the improvement in data quality, tools such as scatter plots, histograms, heatmaps, box plots, and violin plots were employed. The overall process of this study consists of three main stages:

• Initial data analysis without outlier removal (Fig. 1).

• Identification and visualization of outliers, where red triangles represent the outliers, and the symbol “*” denotes the remaining valid data points (Fig. 2).

• Analysis of data quality improvement after outlier removal (Fig. 3).

The Gaussian removal method, or outlier removal using the Gaussian distribution, is a statistical technique used to detect and eliminate outliers in a dataset. This method is based on the assumption that the data follows a normal distribution (Gaussian distribution), and data points that deviate significantly from this distribution are identified as outliers (Fig. 2).

The process of the Gaussian Removal method consists of three main steps. First, the mean \((\mu )\) and standard deviation \((\sigma )\) of the dataset are calculated to identify outliers. Then, data points are identified as outliers if their distance from the mean is greater than one or more times the standard deviation. This distance is usually defined as \(\mu \pm k\alpha\), where \(k\) is a constant value that determines the sensitivity to outliers. Finally, data points that fall outside this range are identified as outliers and removed.

Formulas for calculating mean and standard deviation:

where \({x}_{i}\) are the data points, \(n\) is the number of data points, and \(\mu\) is the mean of the data.

Gaussian removal equation:

Data points identified as outliers follow the following equation:

where \(P\left(x\right)\) is the probability density for data point \(x\), \(\mu\) is the mean of the data, \(\sigma\) is the standard deviation of the data and \(x\) are the data points.

Data points with very low probability values (i.e., significantly distant from the mean) are identified as outliers and removed.

A z-score threshold of ± 3.0 was selected for identifying outliers, corresponding to a 99.7% confidence interval under the assumption of a normal distribution. This threshold is commonly used in statistical analysis to balance the removal of anomalous values without discarding significant but valid data. Its selection was informed by prior studies and tailored to the size and variability of the dataset used in this study.

To scientifically demonstrate the improvement in data quality after outlier removal, various data analysis tools were used, including scatter plots, histograms, heatmaps, box plots, and violin plots. The results of these analyses are presented in Fig. 3.

Machine learning methods

To gain deeper insights into the vertical variation of petrophysical properties, Fig. 4 presents the depth-wise trends of all nine input parameters utilized for machine learning model development. By plotting these well-log features—such as porosity, resistivity, gamma ray, and sonic travel time—as functions of depth, this figure enables the identification of geological transitions, lithological changes, and anomalous zones. Understanding these depth-dependent behaviors is crucial for both feature engineering and improving model interpretability, particularly in heterogeneous formations. These trends also highlight zones of consistent property behavior, which can guide spatially-aware data partitioning and ensure robust training and testing procedures. Depth-wise distribution of the nine well-log input parameters used in the ML models (DEPT, HTNP, RT, CGR, SGR, HCAL, DTCO, RHOZ, and TEMP). These plots reveal the variation of petrophysical properties with depth, highlighting zones of heterogeneity, lithological changes, and potential data anomalies. This visualization aids in understanding reservoir behavior and in developing more geologically consistent and interpretable ML models.

The selection of the five machine learning algorithms—LR, SVM, RF, LSBoost, and Bayesian modeling—was based on their proven effectiveness in regression-based reservoir property prediction and their ability to capture both linear and nonlinear relationships with varying degrees of interpretability and computational efficiency. These models were chosen to cover a broad range of algorithmic families, including parametric (LR, Bayesian), ensemble (RF, LSBoost), and kernel-based approaches (SVM), allowing for a robust comparative analysis across different learning paradigms. While neural networks and deep learning methods have shown great potential in subsurface modeling, they typically require significantly larger datasets, extensive computational resources, and longer training times. Given the structure and volume of our dataset, as well as the objective of developing interpretable and computationally efficient models, we prioritized algorithms that balance performance with clarity and feasibility. Nonetheless, future work may explore deep learning models to further evaluate their predictive capacity in Sw estimation.

To mitigate the risk of data leakage arising from spatial or depth-wise correlation in well-log datasets, the data splitting strategy was revised to incorporate spatial separation. Data points were grouped based on well identifiers and depth intervals, ensuring that samples from the same well or adjacent depths were not simultaneously included in both training and testing sets. This spatially-aware partitioning strategy preserves the geological continuity of the data and provides a more reliable assessment of model generalization. Grouped cross-validation was employed to evaluate model performance under these constraints, thereby improving the robustness of the results.

Bayesian

The Bayesian approach is a fundamental technique in ML that leverages the principles of Bayesian probability to model and learn from data. This method integrates prior knowledge with evidence obtained from the data. At its core lies Bayes’ theorem, which defines the relationship between conditional probabilities (see Fig. 5). Bayes’ Theorem can be written as:

In this equation, \(P\left( {H|D} \right)\) represents the posterior probability, which indicates the likelihood of hypothesis \(H\) given the evidence \(D\). \(P\left( {D|H} \right)\) is the likelihood, showing how probable the data \(D\) is under the assumption that hypothesis \(H\) is true. \(P\left( H \right)\) is the prior probability, representing initial knowledge about hypothesis \(H\), and \(P\left( D \right)\) is the marginal likelihood, acting as a normalizing factor.

In machine learning, \(H\) typically denotes the model or its parameters, while \(D\) represents the training data. The goal is to estimate the posterior probability \(P\left( {D|H} \right)\) to learn the model or its parameters. Bayesian methods are broadly categorized into parametric Bayesian learning and non-parametric Bayesian learning. In parametric learning, the model parameters are assumed to be fixed but unknown. For instance, if \(\theta\) represents the model parameters, the posterior distribution is expressed as:

Conversely, non-parametric Bayesian learning is utilized when the model structure or the number of parameters is not predefined. This method is often applied in models such as Gaussian processes or Bayesian clustering, where the complexity of the model dynamically adapts to fit the data.

One of the key applications of Bayesian methods in ML is prediction. Predictions are made using the posterior expectation. For example, the prediction of \({y}^{*}\) for a new data point \({x}^{*}\) is calculated as:

Here, \(P\left({y}^{*}|{x}^{*},\theta \right)\) is the predictive probability of the output \({y}^{*}\) given the input \({x}^{*}\) and parameters \(\theta\), and \(P\left(D|\theta \right)\) represents the posterior distribution of the parameters. This integral is often computed using methods such as Monte Carlo sampling (MCMC) or other approximation techniques.

Bayesian methods are applied in various models, such as the Naive Bayes Classifier and Bayesian Networks. In the Naive Bayes Classifier, it is assumed that the features are independent of each other, and the probability of a class \(C\) given features \({x}_{1}, {x}_{2}, \dots , {x}_{n}\) is calculated as:

On the other hand, Bayesian Networks leverage the causal relationships between variables. In these models, the nodes represent the variables, while the edges indicate the probabilistic dependencies connecting them.

In this study, the Bayesian method was implemented using a Gaussian prior distribution and the Markov Chain Monte Carlo (MCMC) sampling technique. A burn-in period of 100 iterations and 1000 posterior samples were used to ensure convergence. The model was configured for regression, with hyperparameters selected based on prior studies and refined through performance-based sensitivity analysis. The implementation aimed to evaluate the probabilistic estimation behavior of water saturation under uncertainty, while maintaining interpretability and robustness.

The Bayesian approach offers both strengths and limitations. Among its advantages are the ability to integrate prior knowledge with new data, the provision of probabilistic distributions instead of fixed values, and its suitability for situations with limited data. However, its main challenges include high computational complexity and sensitivity to the choice of the prior distribution.

Linear regression

Linear Regression is a fundamental and commonly applied supervised learning algorithm in machine learning. It aims to capture the relationship between independent variables (features) and a dependent variable (target) by fitting a straight-line equation to the observed data. The main goal is to estimate the dependent variable’s value based on the provided independent variables. The general mathematical formulation of LR is represented as:

where, \(y\) is the dependent variable, \(x_{1} ,x_{2} , \ldots , x_{p}\) are the dependent variable, \(\beta_{0}\) is the intercept, \(\beta_{1} ,\beta_{2} , \ldots , \beta_{p}\) are the coefficients, and \(\varepsilon\) is the error term. In matrix form, it can be written as:

Here, \(y, X, \beta , and \epsilon\) are matrices and vectors representing the data, coefficients, and errors. The main goal of LR is to minimize the error between the predicted and observed values, which is typically measured using the Mean Squared Error (MSE). The formula for MSE is:

where, \({y}_{i}\) represents the actual values, and \({\widehat{y}}_{i}\) represents the predicted values.

The parameters of the model \(\left(\beta \right)\) are estimated using the Ordinary Least Squares (OLS) method. The OLS solution is derived as:

Once the model is trained, predictions for new data points are made using the following formula:

Here, \(x\) is the feature vector of the new data point.

The Linear Regression model was employed using the OLS approach with L2 regularization (ridge regression). A regularization coefficient of α = 0.01 and convergence tolerance of 1e-4 were applied to prevent overfitting while preserving model simplicity. This model served as a baseline to evaluate the linear predictability of Sw from the selected input features and to compare against more complex nonlinear models.

Support vector machine

Support Vector Machines (SVMs) are a powerful supervised learning algorithm widely utilized for tasks such as classification, regression, and outlier detection. They perform well in high-dimensional spaces and can manage both linear and nonlinear classification problems. The key goal of an SVM is to determine a hyperplane that best separates data points belonging to different classes in the feature space.

The central idea of SVM is to determine a hyperplane that maximizes the margin between two classes. This margin is defined as the distance between the hyperplane and the nearest data points from each class, known as support vectors. For a training dataset represented as \(\left\{\left({x}_{i}, {y}_{i}\right)\right\}\), where \({x}_{i}\) is the feature vector and \({y}_{i}\in \left\{-1, +1\right\}\) is the class label, the hyperplane is mathematically defined as:

where, \(w\) represents the weight vector, \(x\) is the input feature vector, and \(b\) is the bias term. The hyperplane acts as the decision boundary, while the support vectors are the data points closest to this boundary.

For linearly separable data, SVM aims to find the optimal hyperplane that maximizes the margin between the two classes. The equations defining the margin boundaries are:

This method guarantees that the hyperplane maximizes the separation between the classes while keeping the nearest points (support vectors) at the margin’s boundary.

The SVM model was configured with a RBF kernel, selected for its ability to capture nonlinear relationships. The regularization parameter was set to 10 and the tolerance was fixed at 0.001, as determined through empirical testing and cross-validation. Feature scaling was performed using z-score standardization prior to training to ensure numerical stability and improve convergence. SVM showed superior performance, particularly in handling the nonlinear and high-dimensional nature of the dataset.

Least squares boosting

Least Squares Boosting (LSBoost) is a ML technique that combines boosting with least squares regression to enhance prediction accuracy. Primarily used for regression tasks, it can also be adapted for classification problems. LSBoost improves upon traditional boosting methods by focusing on minimizing the least squares error of the model. The process starts with a simple initial model, typically the mean of the target values. During each boosting iteration, residuals are calculated by finding the difference between the actual values and the model’s current predictions. A new weak learner, usually a decision tree, is then trained to predict these residuals, aiming to minimize the least squares error. The model is updated by incorporating the scaled predictions of this weak learner, with adjustments made via a learning rate. After a set number of iterations or once the model’s performance stabilizes, the final model is formed by combining the predictions from all weak learners.

The initial prediction is computed as:

where, \({y}_{i}\) is the actual target value for the \(ith\) instance and \(N\) is the total number of instances.

At iteration \(m\), the residuals are calculated as:

where \(\hat{y}_{i} \left( m \right)\) represents the prediction for the \(ith\) instance at iteration \(m\).

A weak learner is then trained on these residuals with the goal of minimizing the least squares error, which is represented as:

where \(f_{m} \left( {x_{i} } \right)\) is the prediction of the weak learner for the \(ith\) instance.

The model is updated by adding the predictions of the weak learner, scaled by a learning rate \(\alpha\):

where \(\alpha\) is the learning rate (also known as the shrinkage parameter).

After \(M\) iterations, the final prediction is given by:

LSBoost was implemented using 50 boosting iterations with a learning rate of 0.1. Decision trees with a maximum depth of 3 were used as base learners to maintain model generalization and avoid overfitting. The model’s parameters were selected based on iterative testing and validation accuracy. LSBoost was specifically included to test the performance of additive ensemble learning in capturing complex dependencies within the Sw data.

Random forest

Random Forest (RF) is a popular ML method known for its effectiveness in classification and regression tasks, particularly when dealing with large and complex datasets. Its strength lies in its ability to reduce variance and avoid overfitting, making it a preferred option. RF is an ensemble technique that combines numerous decision trees, each trained independently on a random subset of the data. The overall prediction is derived by aggregating the outcomes from these individual trees (Fig. 6).

The RF algorithm begins by creating random samples from the training data using bootstrap sampling (sampling with replacement). Each decision tree is trained on a different random subset of the data. To ensure diversity among the trees, a random selection of features is made at each node of the trees. During prediction, the results from all trees are aggregated. For classification tasks, the final prediction is based on majority voting, while for regression tasks, the predictions are averaged.

In classification, the final prediction for a new sample \(x\) is computed using the following formula, where \(T_{1} , T_{2} , \ldots , T_{n}\) represent the decision trees and \(C_{1} , C_{2} , \ldots , C_{n}\) are the possible classes:

Here, \(I\) is an indicator function that equals 1 if \(T_{j} \left( x \right)\) equals \(C_{i}\) and 0 otherwise.

For regression tasks, the final prediction for a new sample \(x\) is calculated as the average of the predictions from all the decision trees:

where \(T_{j} \left( x \right)\) is the predicted value from decision tree \(T_{j}\) for the sample \(x\).

The RF model was constructed with 100 decision trees and a maximum tree depth of 10. Bootstrap sampling was enabled, and a minimum of two samples was required to split each node. This configuration was selected after several trials to balance performance and training efficiency. The model was used to evaluate the effectiveness of bagging-based ensemble techniques in reducing variance and improving prediction accuracy across heterogeneous input features.

Table 2 contains the key control parameters for various machine learning algorithms, which were carefully tuned and optimized to ensure the accuracy and reliability of the results. It includes the values and descriptions of relevant parameters for each algorithm, such as SVM, RF, LSBoost, Bayesian Method, and LR. These settings were selected based on preliminary testing, cross-validation, and best practices reported in the literature to optimize the performance of each algorithm.