Machine Learning-Assisted Prediction of Polymer Flooding Performance in Heterogeneous Carbonate Reservoirs of the Pre-Caspian Basin: A Gradient Boosting and Physics-Informed Approach

Full Text

1. Introduction

The mature oil fields of Kazakhstan’s Pre-Caspian Basin, including Tengiz, Karachaganak, and the legacy fields of the Mangistau Platform such as Uzen and Zhetybai, collectively hold proven reserves exceeding 30 billion barrels of oil equivalent [1]. However, water cut in many of these assets has reached 80–92%, signalling that primary and secondary recovery mechanisms are approaching their practical limits [2, 3]. Enhanced oil recovery through chemical methods—and polymer flooding in particular-has emerged as the most technically mature intervention for improving volumetric sweep efficiency in reservoirs characterised by adverse mobility ratios and high permeability heterogeneity [4].

Polymer flooding operates on the principle of increasing the apparent viscosity of the injected water phase, thereby reducing the water-to-oil mobility ratio :

where  and  are the relative permeabilities to water and oil at their respective saturations, and  and  are the viscosities of water and oil. When , viscous fingering develops and sweep efficiency deteriorates [5]. Partially hydrolysed polyacrylamide (HPAM) is the dominant polymer deployed commercially, yet its performance in high-salinity, high-temperature carbonate environments-conditions characteristic of many Kazakhstani reservoirs-is sensitive to a complex, nonlinear combination of parameters that resist accurate prediction by analytical correlations alone [6, 7].

Full-physics compositional reservoir simulation remains the gold standard for polymer flooding prediction, but it is computationally demanding: a single high-resolution sector model of a Tengiz-type carbonate complex may require 48–120 hours of CPU time per realisation, rendering probabilistic uncertainty analyses and real-time operational decisions practically infeasible [8]. This gap creates a compelling motivation for surrogate modelling approaches that can achieve near-simulation accuracy at a fraction of the computational cost.

Machine learning methods have demonstrated growing utility in EOR prediction tasks over the past decade. Early applications used artificial neural networks to correlate recovery factor with reservoir quality indices [9, 10]. More recently, ensemble tree methods-gradient boosting in particular-have shown superior predictive performance on tabular petroleum engineering datasets due to their native handling of mixed feature types, robustness to outliers, and interpretability through feature importance and SHAP values [11, 12]. Nevertheless, a critical weakness of purely data-driven models is their potential to produce physically implausible outputs: for example, predicting increased injectivity as salinity rises, contrary to the well-established HPAM adsorption-degradation mechanism [13]. Physics-informed ML (PIML) addresses this by embedding physical constraints-either as penalty terms in the loss function, as hard constraints on the output layer, or as physics-derived feature transformations-thereby combining the flexibility of statistical learning with the consistency guarantees of first-principles models [14, 15].

The present study makes the following specific contributions to the literature:

1. A curated, quality-controlled dataset of 412 polymer flooding operations-the largest yet assembled for Kazakhstani and Central Asian carbonate reservoirs-is constructed and made partially available.
2. A physics-informed gradient boosting framework is developed in which a penalty term derived from the Darcy-scale polymer transport equation is integrated into the XGBoost loss function.
3. The model is validated against full compositional simulation using the Uzen carbonate pilot as a case study, demonstrating agreement within ±1.8% absolute in incremental recovery factor.
4. SHAP-based sensitivity analysis provides quantitative ranking of feature importance, yielding physically interpretable and actionable guidance for field engineers.

2. Materials and Methods

2.1 Dataset Construction and Feature Engineering

The primary dataset comprises 412 observations drawn from three sources: (i) 187 records from peer-reviewed SPE and EAGE publications on polymer flooding pilots in carbonate and mixed carbonate-clastic systems worldwide with reservoir conditions comparable to Pre-Caspian Basin formations (temperature 60–130 °C, salinity 30–250 g/L TDS); (ii) 143 records from internal operational reports of TCO, KPO, and OMG polymer injection pilots conducted between 2008 and 2024; and (iii) 82 records from the publically accessible United States Department of Energy EOR database after filtering for chemically comparable systems.

For each observation, 18 input features were extracted (Table 1). The target variables are: (a) incremental oil recovery factor  (%), defined as the difference between the final recovery factor achieved with polymer flooding and the estimated recovery under continued waterflood; and (b) injectivity ratio , the ratio of polymer injection rate to baseline water injection rate at equivalent bottomhole pressure.

Table 1. Input feature set and descriptive statistics (N = 412)

Missing values (6.3% of cells overall) were imputed using an iterative random-forest imputer (IterativeImputer, scikit-learn v1.4.2). The categorical polymer type variable was one-hot encoded. All continuous features were standardised to zero mean and unit variance prior to training of the ANN and SVR baselines; GBR is invariant to monotonic feature scaling but the standardised set was also used for interpretable SHAP value comparison.

2.2 Theoretical Background: Polymer Transport and the Physics-Informed Penalty

The transport of polymer in a porous medium is governed by an advection-dispersion-reaction equation that, in one-dimensional form and neglecting dispersion, reduces to [16]:

where  is the in-situ polymer concentration (kg/m³),  is the adsorbed polymer concentration (kg/kg rock),  is the rock grain density (kg/m³),  is the Darcy flux (m/s), and  is the spatial coordinate. The adsorption isotherm for HPAM on carbonate surfaces is commonly described by the Langmuir relationship:

where  is the maximum adsorption capacity (μg/g rock) and  is the Langmuir coefficient (m³/kg). The apparent viscosity of the polymer solution is modelled by the Carreau-Yasuda equation accounting for shear-thinning:

where  is the zero-shear viscosity (a function of polymer concentration),  is the infinite-shear viscosity,  is the relaxation time,  is the shear rate, and ,  are fitting parameters.

From these governing equations, two physics-derived inequality constraints are formulated for embedding in the ML loss function:

Constraint 1 (Salinity-Injectivity Monotonicity): An increase in divalent ion concentration must not increase polymer injectivity, because divalent cations promote HPAM adsorption and precipitation. Formally: .

Constraint 2 (Concentration-Viscosity Monotonicity): Increasing polymer concentration at fixed molecular weight must not decrease effective viscosity. Formally: , which by the mobility ratio argument implies  up to the adsorption saturation limit.

These constraints are incorporated as soft penalties added to the XGBoost objective function during training. The augmented loss  is defined as:

where  is the mini-batch,  and  are penalty weights tuned via a grid search on the validation set, and the partial derivatives are evaluated numerically by finite difference on the model predictions.

2.3 Model Architecture and Training

Three ML models were trained and benchmarked:

Gradient Boosting Regression (GBR / XGBoost): XGBoost v2.0.3 with the physics-informed augmented loss. Key hyperparameters after Bayesian optimisation (Optuna v3.6): number of estimators ; maximum tree depth ; learning rate ; subsample ratio 0.78; column sample by tree 0.72; L2 regularisation .

Artificial Neural Network (ANN): Four-layer fully connected network with architecture [18 → 128 → 64 → 32 → 1], ReLU activations, dropout rate 0.25, trained with Adam optimiser (lr = 0.001) for 500 epochs with early stopping (patience = 30) on mean squared error.

Support Vector Regression (SVR): Radial basis function kernel; regularisation parameter ; epsilon ; optimised by grid search.

All models were trained on an 80% stratified random split (329 samples) and evaluated on the held-out 20% (83 samples). Stratification was performed on quartiles of  to ensure representative target distribution in both sets. Five-fold cross-validation was used for hyperparameter optimisation. Performance metrics reported include , RMSE (%), and mean absolute error MAE (%).

3. Results

3.1 Predictive Performance Comparison

Table 2 summarises the predictive accuracy of all models on the held-out test set for both target variables.

Table 2. Model performance on held-out test set (N = 83)

The proposed physics-informed GBR model achieves the highest accuracy across both targets, reducing RMSE by 6.1% relative to unconstrained GBR for ΔRF. More significantly, the physics penalty reduced the fraction of predictions violating Constraint 1 (negative IR response to divalent salinity) from 14.2% to 3.1% of test observations (a 78% reduction), and violations of Constraint 2 from 9.7% to 2.4%.

Figure 1: Scatter plot of predicted vs. observed incremental recovery factor (ΔRF, %) for the four models on the 83-sample test set. Each panel shows: x-axis = observed ΔRF (0–25%); y-axis = predicted ΔRF (0–25%); identity line in red; 95% prediction interval bands in shaded grey. Points coloured by reservoir permeability variation coefficient  (blue = low heterogeneity, yellow = high heterogeneity). The physics-informed GBR panel shows the tightest point cloud around the identity line, with notably fewer outliers in the high- (high-heterogeneity) region.

3.2 SHAP Feature Importance Analysis

SHAP values were computed for the physics-informed GBR model using TreeExplainer (SHAP v0.45). Figure 2 presents the global feature importance beeswarm plot for ΔRF prediction.

Figure 2: SHAP beeswarm plot for incremental recovery factor (ΔRF) prediction. x-axis: SHAP value (impact on model output, % ΔRF units); each point represents one observation coloured by feature value (blue = low, red = high). Features ranked top to bottom by mean |SHAP| value. The three dominant features are , , and .

The permeability variation coefficient  emerges as the most influential feature (mean |SHAP| = 3.82%), consistent with the fundamental role of reservoir heterogeneity in controlling polymer sweep efficiency [17]. Polymer concentration  ranks second (2.94%), and formation water salinity  (TDS) third (2.61%). The negative SHAP contribution of high  values confirms the known mechanism of HPAM viscosity loss in saline environments and validates the physical interpretability of the model.

3.3 Sensitivity Analysis: Polymer Concentration and Salinity Effects

To further validate model behaviour against physical expectations, partial dependence plots were generated by varying  (200–3000 mg/L) and  (30–250 g/L) while holding all other features at their median values. Table 3 summarises the predicted ΔRF and IR at selected values.

Table 3. Predicted ΔRF and IR as functions of polymer concentration and salinity (all other features at median values)

The table reveals a diminishing returns relationship between  and ΔRF: incremental gain per 500 mg/L increment decreases from 3.6% (500→1000 mg/L at  = 50 g/L) to 0.8% (2000→2500 mg/L), which is consistent with Langmuir adsorption saturation. The monotonic decrease in both ΔRF and IR with increasing salinity confirms the physics-imposed constraints are operative in the trained model.

Figure 3: 3D surface plot of predicted ΔRF (z-axis, %) as a function of polymer concentration  (x-axis, mg/L) and formation water salinity  (y-axis, g/L TDS). All other features held at median values. Surface coloured by ΔRF magnitude (blue = low, red = high). Isocontour lines plotted at ΔRF intervals of 2%.

3.4 Case Study: Uzen Carbonate Pilot Validation

The Uzen field (Mangistau Region) operates a polymer flooding pilot covering 3 injectors and 8 producers in the Lower Triassic carbonate reservoir (J-XIII horizon, average  = 62 mD,  = 0.74,  = 78°C, TDS = 138 g/L). A benchmark comparison was made between the ML surrogate, the full Eclipse 300 compositional simulation, and the actual field performance recorded over 36 months of polymer injection.

Table 4. Uzen carbonate pilot validation: ML surrogate vs. full simulation vs. field data

The ML surrogate reproduces the Eclipse 300 incremental recovery profile within ±1.8% absolute at all evaluated PVI checkpoints, with a final ΔRF of 12.6% vs. the simulation value of 12.8%. The model also correctly captures the field-observed final ΔRF of 12.1% within its prediction uncertainty band (±2.5% at the 95% confidence level). The computational advantage is dramatic: the surrogate requires less than 1 second per evaluation versus approximately 72 hours for the full compositional simulation, a speedup exceeding five orders of magnitude.

Figure 4: Recovery factor increment ΔRF (%) as a function of pore volumes injected (PVI) for the Uzen carbonate pilot. Three curves are plotted: (1) Eclipse 300 full simulation (solid black line); (2) physics-informed GBR surrogate ± 95% confidence interval (blue line with shaded band); (3) actual field data points (red diamonds, six measurement intervals over 36 months). x-axis: PVI (0 to 2.0); y-axis: cumulative ΔRF (0 to 16%); legend in upper left corner.

4. Discussion

The results demonstrate that the physics-informed gradient boosting surrogate provides a practically useful approximation to full-physics compositional simulation for polymer flooding performance prediction in carbonate reservoirs. Several aspects merit detailed discussion.

4.1 Role of Physics Constraints in Improving Physical Consistency

The 78% reduction in Constraint 1 violations (salinity-injectivity monotonicity) attributable to the physics penalty is particularly noteworthy. Without the penalty, unconstrained GBR occasionally predicted positive injectivity responses to increasing divalent ion concentration in the 50–150 mg/L  range, a physically implausible outcome attributable to spurious correlations in the training data where high-salinity reservoirs sometimes coincidentally received higher concentration polymer slugs. The physics penalty effectively suppresses these artefacts by imposing the thermodynamically consistent monotonic relationship as a soft constraint, rather than relying on the training data alone to encode it.

This finding aligns with the broader PIML literature: Raissi et al. [14] and Karniadakis et al. [15] demonstrated that physics constraints provide the most performance gain precisely when training data are sparse or contain conflating correlations-conditions that characterise petroleum engineering datasets assembled from heterogeneous field operations.

4.2 Dominant Role of Permeability Heterogeneity

The SHAP analysis confirms that  (Dykstra-Parsons coefficient) is the single most influential predictor of ΔRF, consistent with the theoretical framework of Stiles [18] and the extensive simulation study of Sorbie [5]. In highly heterogeneous reservoirs (), the polymer viscosity improvement alone may be insufficient to correct the unequal layer flux distribution; in such cases, cross-linked polymer or polymer-surfactant combinations may be necessary [19]. The ML model implicitly captures this threshold behaviour through the gradient boosting interaction terms, as evidenced by the strongly nonlinear SHAP dependence plot for  (not shown in full here for brevity).

4.3 Economic Screening Implications for Kazakhstani Fields

The 5-orders-of-magnitude speedup over full simulation (sub-second vs. 72 hours) enables uncertainty quantification workflows previously inaccessible for polymer EOR projects in Kazakhstan. For example, a Monte Carlo analysis sampling 10,000 reservoir realisations-routine in modern probabilistic reserves assessment-is completed in less than 3 hours with the surrogate, versus an estimated 72,000 CPU-hours with Eclipse 300. This capability aligns with Kazakhstan’s national EOR development strategy targeting substantial additional oil recovery through chemical flooding over the coming decade.

The model also provides actionable guidance for polymer selection and slug design. For the median Kazakhstani carbonate reservoir conditions in the dataset ( = 87 mD,  = 0.68,  = 82°C, TDS = 112 g/L), the surrogate predicts an optimal polymer concentration window of 1200–1600 mg/L, with diminishing returns above 1800 mg/L. This range is consistent with the laboratory coreflooding results of Seitenova et al. [20] for Tengiz-type dolomite cores at comparable conditions.

4.4 Limitations

Several limitations of the present study should be acknowledged. First, the dataset, while the largest assembled for this class of operations in Central Asian carbonates, contains 25% records from global analogue fields rather than Kazakhstani assets; transfer learning penalties may introduce bias for extreme reservoir conditions outside the Kazakhstani range. Second, the physics constraints embedded are one-dimensional and steady-state; they do not capture transient polymer propagation front dynamics or near-wellbore rheological effects, which may be significant in radial flow geometries. Third, the model was trained and validated at the field-segment scale; extrapolation to individual well predictions requires additional well-level training data currently not available in the dataset. Future work will address these limitations through (i) active learning data acquisition from ongoing Kazakhstan pilot programmes, (ii) incorporation of graph neural network architectures to capture spatial well connectivity, and (iii) extension of physics constraints to capture two-dimensional sweep patterns using proxy streamline models.

5. Conclusions

This study presented and validated a physics-informed gradient boosting surrogate model for predicting polymer flooding performance in heterogeneous carbonate reservoirs of Kazakhstan’s Pre-Caspian Basin. The following conclusions are drawn:

1. The physics-informed GBR model achieved and RMSE = 2.31% for incremental recovery factor prediction on the held-out test set, outperforming ANN, SVR, and unconstrained GBR baselines across all evaluated metrics.
2. Embedding physics-derived inequality constraints (salinity-injectivity monotonicity and concentration-viscosity monotonicity) as soft penalty terms in the XGBoost loss function reduced physically inconsistent predictions by 78% compared to unconstrained GBR, substantially improving model trustworthiness for operational deployment.
3. SHAP analysis identified the Dykstra-Parsons permeability variation coefficient , polymer concentration , and formation water salinity as the three dominant predictors of ΔRF, consistent with established polymer flooding theory.
4. Validation against Eclipse 300 full compositional simulation and actual field measurements on the Uzen J-XIII carbonate pilot confirmed surrogate accuracy within ±1.8% absolute in incremental recovery factor, while delivering a computational speedup exceeding five orders of magnitude.
5. The surrogate framework enables Monte Carlo uncertainty quantification and EOR candidate screening at field portfolio scale within practical computational budgets, directly supporting Kazakhstan’s national EOR development objectives.

About the authors

Adil Almasuly Zhetpis

Email: adiljetpeace@gmail.com
ORCID iD: 0009-0000-0390-5436
Kazakhstan

Vladimir Vitalievich Ryndin

Email: ryndin.v@teachers.tou.edu.kz
ORCID iD: 0000-0002-4248-9516

Mensulu Manarbekovna Kenzhanova

Email: mena-1997@list.ru
ORCID iD: 0009-0004-4696-6908

Inara Kassymkhanovna Omarbekova

Author for correspondence.
Email: inara_19@mail.ru
ORCID iD: 0009-0004-4292-6765

Dinara Gabdulgazisovna Aigozhina

Email: dina7041976@mail.ru
ORCID iD: 0000-0002-5712-5280

Gulnara Gosmanovna Abdullina

Email: gulnara_1277@mail.ru
ORCID iD: 0000-0003-4493-1715

References

Supplementary files