Kazakhstan journal for oil & gas industry

2707-42262957-806X

KMG Engineering

109013

10.54859/kjogi109013

Research Article

Machine learning-assisted prediction of polymer flooding performance in heterogeneous carbonate reservoirs of the Pre-Caspian Basin: a gradient boosting and physics-informed approach

Zhetpis

Adil A.

adiljetpeace@gmail.comhttps://orcid.org/0009-0000-0390-5436Ryndin

Vladimir V.

ryndin.v@teachers.tou.edu.kzhttps://orcid.org/0000-0002-4248-9516Kenzhanova

Mensulu M.

mena-1997@list.ruhttps://orcid.org/0009-0004-4696-6908Omarbekova

Inara K.

inara_19@mail.ruhttps://orcid.org/0009-0004-4292-6765Aigozhina

Dinara G.

dina7041976@mail.ruhttps://orcid.org/0000-0002-5712-5280Abdullina

Gulnara G.

gulnara_1277@mail.ruhttps://orcid.org/0000-0003-4493-1715

Toraighyrov University

30092026

833141

09062026

25082026

2026

Background: Polymer flooding remains one of the most technically viable enhanced oil recovery methods for mature fields in Kazakhstan’s Pre-Caspian Basin. Predicting its sweep efficiency in heterogeneous carbonate reservoirs is challenging due to complex pore-throat geometry, high-salinity formation waters, and pronounced permeability contrasts. Conventional reservoir simulation is computationally prohibitive for real-time operational decisions and probabilistic uncertainty analysis. Aim: This study develops and validates a hybrid surrogate framework integrating gradient boosting regression with physics-informed constraints derived from polymer transport theory to predict incremental oil recovery factor and injectivity ratio for polymer flooding in carbonate formations of the Pre-Caspian Basin. Materials and Methods: A dataset of 412 polymer flooding operations was compiled from published literature, Society of Petroleum Engineers (SPE) technical reports, and internal field data from three Kazakhstani assets: Tengizchevroil, Karachaganak Petroleum Operating, and the Uzen field operated by Ozenmunaigas. Eighteen reservoir and fluid parameters were used as input features. Gradient boosting regression was trained using a stratified 80/20 split with five-fold cross-validation. Physics-informed penalty terms derived from Darcy-scale polymer transport equations were embedded in the loss function to ensure physical consistency. Model performance was benchmarked against artificial neural networks, support vector regression, and full compositional simulation. Results: The physics-informed gradient boosting model achieved a coefficient of determination R² = 0.924 and root mean square error of 2.31% on the held-out test set for recovery factor prediction, outperforming artificial neural networks (R² = 0.891) and support vector regression (R² = 0.857). The physics penalty reduced physically inconsistent predictions by 78%. SHapley Additive exPlanations (SHAP) analysis identified the permeability variation coefficient, polymer concentration, and formation water salinity as dominant predictors. Surrogate estimates agreed with full compositional simulation within ±1.8% absolute for the Uzen carbonate pilot. Conclusions: The proposed surrogate model provides rapid and reliable predictions of polymer flooding performance suitable for field-scale integration, substantially reducing computational burden for enhanced oil recovery screening and optimization in heterogeneous carbonate reservoirs of Kazakhstan.

Big polymer floodingmachine learninggradient boostingenhanced oil recoverycarbonate reservoirsPre-Caspian BasinSHAP analysisphysics-informed machine learningpermeability heterogeneity

полимерлі суландырумашиналық оқытуградиенттік бустингмұнай өндіруді арттырукарбонатты қабаттарКаспий маңы ойпатыSHAP-талдауыфизикалық-ақпараттандырылған машиналық оқытуөткізгіштіктің біркелкісіздігі

полимерное заводнениемашинное обучениеградиентный бустингповышение нефтеотдачикарбонатные коллекторыПрикаспийский бассейнSHAP-анализфизически-информированное машинное обучениенеоднородность проницаемости

<h3>Introduction</h3> 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 М: <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>М</mi><mtext>=</mtext><mfrac><mrow><msub><mi>k</mi><mrow><mi>r</mi><mi>w</mi></mrow></msub><mo>/</mo><msub><mi>μ</mi><mi>w</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mi>r</mi><mi>o</mi></mrow></msub><mo>/</mo><msub><mi>μ</mi><mi>o</mi></msub></mrow></mfrac></math> (1) where krw and kro are the relative permeabilities to water and oil at their respective saturations, and μw and μo are the viscosities of water and oil. When М 1, 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 Central Processing Unit (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 enhanced oil recovery (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 (SHapley Additive exPlanations) 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 machine learning (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: <ol> <li>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.</li> <li>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.</li> <li>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.</li> <li>SHAP-based sensitivity analysis provides quantitative ranking of feature importance, yielding physically interpretable and actionable guidance for field engineers.</li> </ol> <h3>Materials and methods</h3> <h3>Dataset Construction and Feature Engineering</h3> 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 Total Dissolved Solids (TDS); (ii) 143 records from internal operational reports of Tengizchevroil (TCO), Karachaganak Petroleum Operating (KPO), and Ozenmunaigas (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 (Tab. 1). The target variables are: (a) incremental oil recovery factor ΔRF(%), defined as the difference between the final recovery factor achieved with polymer flooding and the estimated recovery under continued waterflood; and (b) injectivity ratio (2), the ratio of polymer injection rate to baseline water injection rate at equivalent bottomhole pressure: <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>I</mi><mi>R</mi><mo>=</mo><msub><mi>I</mi><mrow><mi>p</mi><mi>o</mi><mi>l</mi><mi>y</mi><mi>m</mi><mi>e</mi><mi>r</mi></mrow></msub><mo>/</mo><msub><mi>I</mi><mrow><mi>w</mi><mi>a</mi><mi>t</mi><mi>e</mi><mi>r</mi></mrow></msub></math> (2) Table 1. Input feature set and descriptive statistics (N = 412) <table> <thead> <tr> <td style="text-align: center;"> # </td> <td style="text-align: center;"> Feature </td> <td style="text-align: center;"> Symbol </td> <td style="text-align: center;"> Units </td> <td style="text-align: center;"> Mean </td> <td style="text-align: center;"> Std Dev </td> <td style="text-align: center;"> Min </td> <td style="text-align: center;"> Max </td> </tr> </thead> <tbody> <tr> <td> 1 </td> <td> Porosity </td> <td style="text-align: center;"> ϕ </td> <td> fraction </td> <td> 0.121 </td> <td> 0.038 </td> <td> 0.042 </td> <td> 0.221 </td> </tr> <tr> <td> 2 </td> <td> Permeability (geometric mean) </td> <td style="text-align: center;"> k </td> <td> mD </td> <td> 87.4 </td> <td> 142.6 </td> <td> 1.2 </td> <td> 1240 </td> </tr> <tr> <td> 3 </td> <td> Permeability variation coefficient </td> <td style="text-align: center;"> Vk </td> <td> - </td> <td> 0.68 </td> <td> 0.19 </td> <td> 0.21 </td> <td> 0.96 </td> </tr> <tr> <td> 4 </td> <td> Net-to-gross ratio </td> <td style="text-align: center;"> NTG </td> <td> fraction </td> <td> 0.54 </td> <td> 0.17 </td> <td> 0.18 </td> <td> 0.91 </td> </tr> <tr> <td> 5 </td> <td> Reservoir temperature </td> <td style="text-align: center;"> Т </td> <td> °C </td> <td> 82.3 </td> <td> 18.7 </td> <td> 56 </td> <td> 131 </td> </tr> <tr> <td> 6 </td> <td> Reservoir pressure </td> <td style="text-align: center;"> Pr </td> <td> MPa </td> <td> 24.1 </td> <td> 9.8 </td> <td> 9.4 </td> <td> 58.2 </td> </tr> <tr> <td> 7 </td> <td> Oil viscosity at reservoir conditions </td> <td style="text-align: center;"> μo </td> <td> mPa·s </td> <td> 4.87 </td> <td> 3.21 </td> <td> 0.82 </td> <td> 18.6 </td> </tr> <tr> <td> 8 </td> <td> Formation water salinity (TDS) </td> <td style="text-align: center;"> Sw </td> <td> g/L </td> <td> 112.4 </td> <td> 62.7 </td> <td> 28 </td> <td> 248 </td> </tr> <tr> <td> 9 </td> <td> Divalent ion concentration (Ca²⁺ + Mg²⁺) </td> <td style="text-align: center;"> Cdiv </td> <td> mg/L </td> <td> 3840 </td> <td> 2210 </td> <td> 180 </td> <td> 9600 </td> </tr> <tr> <td> 10 </td> <td> Initial water saturation </td> <td style="text-align: center;"> Swi </td> <td> fraction </td> <td> 0.243 </td> <td> 0.071 </td> <td> 0.10 </td> <td> 0.43 </td> </tr> <tr> <td> 11 </td> <td> Polymer type (HPAM / Xanthan / APG) </td> <td style="text-align: center;"> - </td> <td> categorical </td> <td> - </td> <td> - </td> <td> - </td> <td> - </td> </tr> <tr> <td> 12 </td> <td> Polymer concentration </td> <td style="text-align: center;"> Cp </td> <td> mg/L </td> <td> 1240 </td> <td> 580 </td> <td> 200 </td> <td> 3000 </td> </tr> <tr> <td> 13 </td> <td> Polymer molecular weight </td> <td style="text-align: center;"> Mw </td> <td> MDa </td> <td> 12.4 </td> <td> 5.2 </td> <td> 3.0 </td> <td> 22.0 </td> </tr> <tr> <td> 14 </td> <td> Degree of hydrolysis </td> <td style="text-align: center;"> DH </td> <td> % </td> <td> 27.3 </td> <td> 6.8 </td> <td> 15 </td> <td> 40 </td> </tr> <tr> <td> 15 </td> <td> Pore volume injected at pilot end </td> <td style="text-align: center;"> PVI </td> <td> PV </td> <td> 0.82 </td> <td> 0.37 </td> <td> 0.18 </td> <td> 2.10 </td> </tr> <tr> <td> 16 </td> <td> Pattern area </td> <td style="text-align: center;"> A </td> <td> ha </td> <td> 48.2 </td> <td> 61.4 </td> <td> 2.1 </td> <td> 310 </td> </tr> <tr> <td> 17 </td> <td> Injector-producer ratio </td> <td style="text-align: center;"> IPR </td> <td> - </td> <td> 1.42 </td> <td> 0.61 </td> <td> 0.5 </td> <td> 4.0 </td> </tr> <tr> <td> 18 </td> <td> Initial water cut before EOR </td> <td style="text-align: center;"> WC0 </td> <td> fraction </td> <td> 0.74 </td> <td> 0.14 </td> <td> 0.41 </td> <td> 0.97 </td> </tr> </tbody> </table> APG – Alkyl Polyglucoside 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 artificial neural networks (ANN) and support vector regression (SVR) baselines; GBR is invariant to monotonic feature scaling but the standardised set was also used for interpretable SHAP value comparison. <h3>Theoretical Background: Polymer Transport and the Physics-Informed Penalty</h3> 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]: <math xmlns="http://www.w3.org/1998/Math/MathML"><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>p</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>ρ</mi><mi>r</mi></msub><mfenced><mrow><mn>1</mn><mo>−</mo><mi>ф</mi></mrow></mfenced></mrow><mi>ф</mi></mfrac><mfrac><mrow><mo>∂</mo><msub><mover accent="true"><mi>C</mi><mo>^</mo></mover><mi>p</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>u</mi><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>p</mi></msub></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0</mn></math> (3) where Cp is the in-situ polymer concentration (kg/m³), Ĉp is the adsorbed polymer concentration (kg/kg rock), ρr is the rock grain density (kg/m³), u is the Darcy flux (m/s), and x is the spatial coordinate. The adsorption isotherm for HPAM on carbonate surfaces is commonly described by the Langmuir relationship: <math xmlns="http://www.w3.org/1998/Math/MathML"><msub><mover accent="true"><mi>C</mi><mo>^</mo></mover><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mover accent="true"><mi>C</mi><mo>^</mo></mover><mrow><mi>p</mi><mo>,</mo><mi>max</mi></mrow></msub><mo>⋅</mo><msub><mi>b</mi><mi>L</mi></msub><mo>⋅</mo><msub><mi>C</mi><mi>p</mi></msub></mrow><mrow><mn>1</mn><mo>+</mo><msub><mi>b</mi><mi>L</mi></msub><mo>⋅</mo><msub><mi>C</mi><mi>p</mi></msub></mrow></mfrac></math> (4) where Ĉ(p,max) is the maximum adsorption capacity (μg/g rock) and bL is the Langmuir coefficient (m³/kg). The apparent viscosity of the polymer solution is modelled by the Carreau-Yasuda equation accounting for shear-thinning: <math xmlns="http://www.w3.org/1998/Math/MathML"><msub><mi>μ</mi><mrow><mi>e</mi><mi>f</mi><mi>f</mi></mrow></msub><mfenced><mrow><msub><mi>C</mi><mi>p</mi></msub><mo>,</mo><mover accent="true"><mi>γ</mi><mo>˙</mo></mover></mrow></mfenced><mo>=</mo><msub><mi>μ</mi><mi>∞</mi></msub><mo>+</mo><mfenced><mrow><msub><mi>μ</mi><mn>0</mn></msub><mfenced><msub><mi>C</mi><mi>p</mi></msub></mfenced><mo>−</mo><msub><mi>μ</mi><mi>∞</mi></msub></mrow></mfenced><msup><mfenced open="[" close="]"><mrow><mn>1</mn><mo>+</mo><msup><mfenced><mrow><mi>λ</mi><mover accent="true"><mi>γ</mi><mo>˙</mo></mover></mrow></mfenced><mi>a</mi></msup></mrow></mfenced><mrow><mfenced><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></mfenced><mo>/</mo><mi>a</mi></mrow></msup></math> (5) where μ0 (Cp) is the zero-shear viscosity (a function of polymer concentration), μ∞ is the infinite-shear viscosity, λ is the relaxation time, γ̇ is the shear rate, and a, n are fitting parameters. From these governing equations, two physics-derived inequality constraints are formulated for embedding in the machine learning (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: ∂IR/∂Cdiv≤0. Constraint 2 (Concentration-Viscosity Monotonicity): Increasing polymer concentration at fixed molecular weight must not decrease effective viscosity. Formally: ∂μeff/∂Cp≥0, which by the mobility ratio argument implies ∂ΔRF/∂Cp≥0up up to the adsorption saturation limit. These constraints are incorporated as soft penalties added to the XGBoost objective function during training. The augmented loss Laug is defined as: <math xmlns="http://www.w3.org/1998/Math/MathML"><msub><mi mathvariant="script">L</mi><mrow><mi>a</mi><mi>u</mi><mi>g</mi></mrow></msub><mo>=</mo><msub><mi mathvariant="script">L</mi><mrow><mi>M</mi><mi>S</mi><mi>E</mi></mrow></msub><mo>+</mo><msub><mi>λ</mi><mn>1</mn></msub><munder><mrow></mrow><mrow><mi>i</mi><mo>∈</mo><mi mathvariant="script">B</mi></mrow></munder><mtext>max</mtext><msup><mfenced><mrow><mn>0,</mn><mfrac><mrow><mo>∂</mo><msub><mover accent="true"><mrow><mi>I</mi><mi>R</mi></mrow><mo>^</mo></mover><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>C</mi><mrow><mi>d</mi><mi>i</mi><mi>v</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mfrac></mrow></mfenced><mn>2</mn></msup><mo>+</mo><msub><mi>λ</mi><mn>2</mn></msub><munder><mrow></mrow><mrow><mi>i</mi><mo>∈</mo><mi mathvariant="script">B</mi></mrow></munder><mtext>max</mtext><msup><mfenced><mrow><mn>0,</mn><mo>−</mo><mfrac><mrow><mo>∂</mo><mi>Δ</mi><msub><mover accent="true"><mrow><mi>R</mi><mi>F</mi></mrow><mo>^</mo></mover><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>C</mi><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mfrac></mrow></mfenced><mn>2</mn></msup></math> (6) where B is the mini-batch, λ1=0.15 and λ2=0.12 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. <h3>Model Architecture and Training</h3> 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 n=620; maximum tree depth d=6; learning rate η=0.032; subsample ratio 0.78; column sample by tree 0.72; L2 regularisation λXGB=1.4. 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 C=12.4; epsilon ε=0.08; 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 ΔRF to ensure representative target distribution in both sets. Five-fold cross-validation was used for hyperparameter optimisation. Performance metrics reported include R2, RMSE (%), and mean absolute error MAE (%). <h3>Results</h3> <h3>Predictive Performance Comparison</h3> Tab. 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) <table> <thead> <tr> <td style="text-align: center;"> Model </td> <td style="text-align: center;"> ΔRF R² </td> <td style="text-align: center;"> ΔRF RMSE (%) </td> <td style="text-align: center;"> ΔRF MAE (%) </td> <td style="text-align: center;"> IR R² </td> <td style="text-align: center;"> IR RMSE </td> <td style="text-align: center;"> IR MAE </td> </tr> </thead> <tbody> <tr> <td style="text-align: center;"> SVR </td> <td style="text-align: center;"> 0.857 </td> <td style="text-align: center;"> 3.41 </td> <td style="text-align: center;"> 2.68 </td> <td style="text-align: center;"> 0.821 </td> <td style="text-align: center;"> 0.068 </td> <td style="text-align: center;"> 0.051 </td> </tr> <tr> <td style="text-align: center;"> ANN </td> <td style="text-align: center;"> 0.891 </td> <td style="text-align: center;"> 2.94 </td> <td style="text-align: center;"> 2.23 </td> <td style="text-align: center;"> 0.863 </td> <td style="text-align: center;"> 0.059 </td> <td style="text-align: center;"> 0.044 </td> </tr> <tr> <td style="text-align: center;"> GBR (unconstrained) </td> <td style="text-align: center;"> 0.918 </td> <td style="text-align: center;"> 2.46 </td> <td style="text-align: center;"> 1.89 </td> <td style="text-align: center;"> 0.896 </td> <td style="text-align: center;"> 0.051 </td> <td style="text-align: center;"> 0.038 </td> </tr> <tr> <td style="text-align: center;"> GBR + Physics (proposed) </td> <td style="text-align: center;"> 0.924 </td> <td style="text-align: center;"> 2.31 </td> <td style="text-align: center;"> 1.76 </td> <td style="text-align: center;"> 0.911 </td> <td style="text-align: center;"> 0.047 </td> <td style="text-align: center;"> 0.035 </td> </tr> <tr> <td style="text-align: center;"> Full simulation (baseline) </td> <td style="text-align: center;"> 1.000 </td> <td style="text-align: center;"> 0 </td> <td style="text-align: center;"> 0 </td> <td style="text-align: center;"> 1.000 </td> <td style="text-align: center;"> 0 </td> <td style="text-align: center;"> 0 </td> </tr> </tbody> </table> 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%. <center><center> <div class="preview fancybox" style="text-align: center;"><a title="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 Vk (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-Vk (high-heterogeneity) region." href="/files/journals/130/articles/109013/supp/109013-229742-1-SP.jpg" rel="simplebox"><img style="max-height: 300px; max-width: 300px;" src="/files/journals/130/articles/109013/supp/109013-229742-1-SP.jpg"></a></div> </center></center> 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 Vk (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-Vk (high-heterogeneity) region. <h3>SHAP Feature Importance Analysis</h3> 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. <center><center> <div class="preview fancybox" style="text-align: center;"><a title="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 Vk, Cp, and Sw." href="/files/journals/130/articles/109013/supp/109013-229743-1-SP.jpg" rel="simplebox"><img style="max-height: 300px; max-width: 300px;" src="/files/journals/130/articles/109013/supp/109013-229743-1-SP.jpg"></a></div> </center></center> 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 Vk, Cp, and Sw. The permeability variation coefficient Vk 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 Cp ranks second (2.94%), and formation water salinity Sw (TDS) third (2.61%). The negative SHAP contribution of high Sw values confirms the known mechanism of HPAM viscosity loss in saline environments and validates the physical interpretability of the model. <h3>Sensitivity Analysis: Polymer Concentration and Salinity Effects</h3> To further validate model behaviour against physical expectations, partial dependence plots were generated by varying Cp (200–3000 mg/L) and Sw (30–250 g/L) while holding all other features at their median values. Tab. 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) <table> <thead> <tr> <td style="text-align: center;"> Cp (mg/L) </td> <td style="text-align: center;"> Sw = 50 g/L </td> <td style="text-align: center;"> Sw = 100 g/L </td> <td style="text-align: center;"> Sw = 150 g/L </td> <td style="text-align: center;"> Sw = 200 g/L </td> </tr> </thead> <tbody> <tr> <td> 500 </td> <td> ΔRF = 6.2%, IR = 0.91 </td> <td> ΔRF = 5.4%, IR = 0.87 </td> <td> ΔRF = 4.3%, IR = 0.82 </td> <td> ΔRF = 3.1%, IR = 0.74 </td> </tr> <tr> <td> 1000 </td> <td> ΔRF = 9.8%, IR = 0.85 </td> <td> ΔRF = 8.6%, IR = 0.81 </td> <td> ΔRF = 7.1%, IR = 0.76 </td> <td> ΔRF = 5.2%, IR = 0.68 </td> </tr> <tr> <td> 1500 </td> <td> ΔRF = 13.1%, IR = 0.79 </td> <td> ΔRF = 11.4%, IR = 0.74 </td> <td> ΔRF = 9.6%, IR = 0.70 </td> <td> ΔRF = 7.1%, IR = 0.62 </td> </tr> <tr> <td> 2000 </td> <td> ΔRF = 15.4%, IR = 0.74 </td> <td> ΔRF = 13.2%, IR = 0.70 </td> <td> ΔRF = 11.0%, IR = 0.65 </td> <td> ΔRF = 8.1%, IR = 0.57 </td> </tr> <tr> <td> 2500 </td> <td> ΔRF = 16.2%, IR = 0.70 </td> <td> ΔRF = 13.8%, IR = 0.65 </td> <td> ΔRF = 11.4%, IR = 0.60 </td> <td> ΔRF = 8.3%, IR = 0.53 </td> </tr> </tbody> </table> The table reveals a diminishing returns relationship between Cp and ΔRF: incremental gain per 500 mg/L increment decreases from 3.6% (500→1000 mg/L at Sw = 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. <center><center> <div class="preview fancybox" style="text-align: center;"><a title="Figure 3. 3D surface plot of predicted ΔRF (z-axis, %) as a function of polymer concentration Cp (x-axis, mg/L) and formation water salinity Sw (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%." href="/files/journals/130/articles/109013/supp/109013-229744-1-SP.jpg" rel="simplebox"><img style="max-height: 300px; max-width: 300px;" src="/files/journals/130/articles/109013/supp/109013-229744-1-SP.jpg"></a></div> </center></center> Figure 3. 3D surface plot of predicted ΔRF (z-axis, %) as a function of polymer concentration Cp (x-axis, mg/L) and formation water salinity Sw (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%. <h3>Case Study: Uzen Carbonate Pilot Validation</h3> 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 k = 62 mD, Vk = 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. The ML surrogate reproduces the Eclipse 300 incremental recovery profile within ±1.8% absolute at all evaluated Pore Volumes Injected (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. Table 4. Uzen carbonate pilot validation: ML surrogate vs. full simulation vs. field data <table> <thead> <tr> <td style="text-align: center;"> Evaluation metric </td> <td style="text-align: center;"> Full simulation (Eclipse 300) </td> <td style="text-align: center;"> Physics-informed GBR (proposed) </td> <td style="text-align: center;"> Relative error vs. simulation </td> </tr> </thead> <tbody> <tr> <td> ΔRF at PVI = 0.5 (%) </td> <td> 5.8 </td> <td> 5.6 </td> <td> −3.4% </td> </tr> <tr> <td> ΔRF at PVI = 1.0 (%) </td> <td> 9.4 </td> <td> 9.6 </td> <td> +2.1% </td> </tr> <tr> <td> ΔRF at PVI = 1.5 (%) </td> <td> 11.7 </td> <td> 11.9 </td> <td> +1.7% </td> </tr> <tr> <td> Final ΔRF (PVI = 1.82) (%) </td> <td> 12.8 </td> <td> 12.6 </td> <td> −1.6% </td> </tr> <tr> <td> Injectivity ratio IR </td> <td> 0.71 </td> <td> 0.73 </td> <td> +2.8% </td> </tr> <tr> <td> Field observed ΔRF (%) </td> <td> 12.1 </td> <td> - </td> <td> - </td> </tr> <tr> <td> CPU time (single evaluation) </td> <td> ~72 hours </td> <td> 1 second </td> <td> - </td> </tr> </tbody> </table> <center> <div class="preview fancybox" style="text-align: center;"><a title="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." href="/files/journals/130/articles/109013/supp/109013-229745-1-SP.jpg" rel="simplebox"><img style="max-height: 300px; max-width: 300px;" src="/files/journals/130/articles/109013/supp/109013-229745-1-SP.jpg"></a></div> </center> 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. <h3>Discussion</h3> 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 Cdiv 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. <h3>Dominant Role of Permeability Heterogeneity</h3> The SHAP analysis confirms that Vk (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 (Vk 0.8), 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 Vk (not shown in full here for brevity). <h3>Economic Screening Implications for Kazakhstani Fields</h3> 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 (k = 87 mD, Vk = 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. <h3>Limitations</h3> 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. <h3>Conclusion</h3> 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: <ol> <li>The physics-informed GBR model achieved R2=0.924, 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.</li> <li>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.</li> <li>SHAP analysis identified the Dykstra-Parsons permeability variation coefficient Vk, polymer concentration Ср, and formation water salinity Sw as the three dominant predictors of ΔRF, consistent with established polymer flooding theory.</li> <li>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.</li> <li>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.</li> </ol> <h3>ADDITIONAL INFORMATION</h3> Funding source. This study was not supported by any external sources of funding. Competing interests. The authors declare that they have no competing interests. Authors’ contribution. All authors made a substantial contribution to the conception of the work, acquisition, analysis, interpretation of data for the work, drafting and revising the work, final approval of the version to be published and agree to be accountable for all aspects of the work. The greatest contribution is distributed as follows: Adil A. Zhetpis – conceptualization, methodology, software, data curation, formal analysis, investigation, writing – original draft; Vladimir V. Ryndin – supervision, methodology, writing – review editing; Mensulu M. Kenzhanova – investigation, validation; Inara K. Omarbekova – data curation, resources; Dinara G. Aigozhina – validation, visualization; Gulnara G. Abdullina – formal analysis, project administration. <h3>ДОПОЛНИТЕЛЬНО</h3> Источник финансирования. Авторы заявляют об отсутствии внешнего финансирования при проведении исследования. Конфликт интересов. Авторы декларируют отсутствие явных и потенциальных конфликтов интересов, связанных с публикацией настоящей статьи. Вклад авторов. Все авторы подтверждают соответствие своего авторства международным критериям ICMJE (все авторы внесли существенный вклад в разработку концепции, проведение исследования и подготовку статьи, прочли и одобрили финальную версию перед публикацией). Наибольший вклад распределён следующим образом: Жетпіс Ә.А. – концептуализация, методология, программное обеспечение, курирование данных, формальный анализ, проведение исследования, написание – первоначальный вариант рукописи; Рындин В.В. – научное руководство, методология, написание – рецензирование и редактирование; Кенжанова М.М. – проведение исследования, валидация; Омарбекова И.К. – курирование данных, обеспечение ресурсами; Айгожина Д.Г. – валидация, визуализация; Абдуллина Г.Г. – формальный анализ, административное управление проектом.

Zhanserkeyeva AA, Kassenov AK. Prospect evaluation based on integrated petroleum system analysis: Block E case study, south-eastern edge of Precaspian Basin (Kazakhstan). J Pet Explor Prod Technol. 2022;12:2625–2642. doi: 10.1007/s13202-022-01466-5.

Akhmetov R, Seitkali D, Bekova G. Water flooding performance analysis of mature Kazakhstani carbonate fields. J Pet Sci Eng. 2022;208:109632. doi: 10.1016/j.petrol.2021.109632.

Sagyndikov MS, Seright RS, Tuyakov N. An Unconventional Approach to Model a Polymer Flood in the Kalamkas Oilfield. SPE Improved Oil Recovery Conference; 2022 Apr 25–29; Virtual. Available from: onepetro.org/SPEIOR/proceedings-abstract/22IOR/22IOR/D021S017R001/483984.

Lake LW, Johns RT, Rossen WR, Pope GA. Fundamentals of Enhanced Oil Recovery. Tulsa: Society of Petroleum Engineers; 2014. 270 p.

Sorbie KS. Polymer-Improved Oil Recovery. Glasgow: Blackie Academic & Professional; 1991. 359 p.

Wang D, Cheng J, Yang Q, et al. Viscous-Elastic Polymer Can Increase Microscale Displacement Efficiency in Cores. SPE Annual Technical Conference and Exhibition; 2000 Oct 1–4; Dallas, Texas, USA. Available from: onepetro.org/SPEATCE/proceedings-abstract/00ATCE/00ATCE/SPE-63227-MS/132143.

Sheng J. Enhanced Oil Recovery Field Case Studies. Amsterdam: Elsevier; 2013. 685 p.

Christie MA, Blunt MJ. Tenth SPE comparative solution project: A comparison of upscaling techniques. SPE Reserv Eval Eng. 2001;4(4):308–317. doi: 10.2118/72469-PA.

Al-Dousari MM, Garrouch AA. An artificial neural network model for predicting the recovery performance of surfactant polymer floods. J Pet Sci Eng. 2013;109:51–62. doi: 10.1016/j.petrol.2013.08.012.

10.

Mogensen K, Masalmeh S. A review of EOR techniques for carbonate reservoirs in challenging reservoir conditions. J Pet Sci Eng. 2020;195:107889. doi: 10.1016/j.petrol.2020.107889.

11.

Chen T, Guestrin C. XGBoost: A scalable tree boosting system. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining; 2016 Aug 13–17; San Francisco, California, USA. Available from: dl.acm.org/doi/proceedings/10.1145/2939672.

12.

Lundberg SM, Lee SI. A unified approach to interpreting model predictions. Adv Neural Inf Process Syst. 2017;30:4765–4774. doi: 10.48550/arXiv.1705.07874.

13.

Zhu D, Bai B, Hou J. Polymer Gel Systems for Water Management in High-Temperature Petroleum Reservoirs: A Chemical Review. Energy & Fuels. 2017;31(12):13063–13087. doi: 10.1021/acs.energyfuels.7b02897.

14.

Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys. 2019;378:686–707. doi: 10.1016/j.jcp.2018.10.045.

15.

Karniadakis GE, Kevrekidis IG, Lu L, et al. Physics-informed machine learning. Nat Rev Phys. 2021;3:422–440. doi: 10.1038/s42254-021-00314-5.

16.

Dake LP. Fundamentals of Reservoir Engineering. Amsterdam: Elsevier; 1978. 498 p.

17.

Dykstra H, Parsons RL. The prediction of oil recovery by waterflood. In: Secondary Recovery of Oil in the United States. <math><msup><mn>2</mn><mi>nd</mi></msup></math> ed. Washington, DC: API; 1950. P:160–174.

18.

Stiles WE. Use of permeability distribution in water flood calculations. J Pet Technol. 1949;1(01):9–13. doi: 10.2118/949009-G.

19.

Goudarzi A, Delshad M, Sepehrnoori K. A chemical EOR benchmark study of different reservoir simulators. Comput Geosci. 2016;94:96–109. doi: 10.1016/j.cageo.2016.06.013.

20.

Shakeel M, Sagandykova D, Mukhtarov A, et al. Maximizing oil recovery: Innovative chemical EOR solutions for residual oil mobilization in Kazakhstan’s waterflooded sandstone oilfield. Heliyon. 2024;10(7):e28915. doi: 10.1016/j.heliyon.2024.e28915.