The key challenges in heavy crude oil pipeline transport are associated with the accurate prediction of rheological properties, hydraulic resistance, and flow assurance strategies. Heavy oils typically exhibit pronounced non-Newtonian behavior, including shear-thinning and yield-stress effects, making conventional viscosity correlations inadequate for describing flow under dynamic temperature and compositional conditions (Souas et al., 2020; Falconi et al., 2018; Saleh et al., 2021; Robidas and Gogoi, 2020). Hydraulic resistance is another major concern: the high viscosity of these oils causes significant frictional losses, requiring increased pumping capacity and leading to elevated operational costs (Souas et al., 2020; Jing et al., 2023; Saleh et al., 2021). Moreover, in cold climates, precipitation of waxes and asphaltenes further complicates flow, threatening pipeline integrity and necessitating advanced flow assurance measures (Jing et al., 2023; Lyu et al., 2019).

A variety of viscosity and rheological models have been proposed for crude oil transport, yet most face significant limitations when applied to heavy oils. Arrhenius- and Refutas-type models are widely used for temperature–viscosity relationships but are often unable to account for non-Newtonian effects in oil mixtures (Robidas and Gogoi, 2020; Jing et al., 2020). Bingham and Herschel–Bulkley models are useful for fluids with yield stress but may fail to capture the full complexity of heavy oil rheology under varying dilution and thermal conditions (Falconi et al., 2018; Saleh et al., 2021; Hoshyargar and Ashrafizadeh, 2013). Similarly, Koval and power-law models have been applied to emulsions and multiphase systems, but their predictive power decreases at extreme shear rates and compositions (Saleh et al., 2021; Hoshyargar and Ashrafizadeh, 2013; Pilehvari et al., 1988). Finally, empirical correlations are often developed for light and medium crudes, making them unsuitable for predicting flow behavior of heavy and non-Newtonian oils, particularly under cold conditions (Silva-Oliver et al., 2020; Robidas and Gogoi, 2020; Jing et al., 2020). These shortcomings highlight the need for more comprehensive approaches that integrate temperature, concentration, and molecular interaction effects.

The study of changes in the rheological properties of non-Newtonian crude oils during blending represents a challenging task, since crude oil is a multicomponent fluid with a highly diverse chemical composition.

Additive models traditionally applied to petroleum systems do not always provide sufficiently accurate predictions of property variations in mixtures containing non-Newtonian crude oils. This limitation often leads to increased pumping costs and reduced energy efficiency of pipeline transportation.

Non-Newtonian crude oil is a hydrocarbon fluid whose rheological behavior deviates from the classical Newtonian law, exhibiting effective viscosity dependence not only on temperature T and pressure P , but also on shear rate γ ˙ , deformation time t , and shear stress τ . Thus, in a generalized form, the effective viscosity may be expressed as Equation 1: η = f T , P , γ ˙ , t , τ . (1)

To predict the viscosity of binary crude oil mixtures, a variety of mathematical models have been developed. Table 1 summarizes the most widely used formulations, including both theoretical and empirical approaches. These models can be broadly classified into three groups according to their underlying principles and degree of empirical adjustment.

The Arrhenius, Bingham, Koval, and Kendall–Monroe models are grounded in physical principles that assume ideal or simplified mixing of components. While mathematically straightforward and computationally efficient, they fail to capture the complexity of real systems. In particular, they neglect intermolecular interactions, structural effects, and molecular associations. As a result, these models are inadequate when the component viscosities differ significantly, since they cannot describe the nonlinear effects observed in real crude oil mixtures or account for density and temperature influences.

The Wallace–Henry, Al-Basrah, Chevron, Shan–Peng (No. 1 and No. 2), and Walther models extend theoretical constructs—typically involving logarithmic or mean-value relationships—by incorporating empirical correction terms or employing transformed equations calibrated against experimental data. Compared with purely theoretical models, they offer improved predictive accuracy by considering the influence of fluid density, interaction mechanisms, and correction coefficients. However, their use is generally restricted to specific types of crude oils or mixtures for which they were originally developed.

Finally, the Cragoe, Maxwell, Refutas, Parkash, and Latour/Miadonye models are entirely empirical in nature, derived from large sets of experimental measurements. They are often expressed as complex exponential or double-logarithmic functions obtained through regression analysis. These models can successfully reproduce viscosity behavior across a wide range of crude oils and blending ratios, making them particularly suitable for engineering calculations. Nevertheless, their application requires precise input parameters and a relatively high degree of computational accuracy.

In general, these mathematical models describe mixture viscosity as a function of the viscosities of the constituent components, while neglecting the influence of temperature, pressure, and shear rate variations.

To overcome the limitations of conventional models, recent research has introduced new approaches that integrate both theoretical developments and practical engineering solutions. Temperature- and dilution-dependent viscosity models have improved prediction accuracy for transport limits and phase inversion points in oil mixtures (Silva-Oliver et al., 2020; Jing et al., 2020). Complementary to this, computational fluid dynamics (CFD) and artificial intelligence (AI) techniques are increasingly used to simulate complex non-Newtonian flow behaviors and optimize pipeline design and operation (Saleh et al., 2021; Jing et al., 2020; Al-Maamari et al., 2015).

On the applied side, chemical and mechanical innovations have demonstrated considerable potential to enhance flow assurance and reduce energy consumption. The development of biodegradable additives, advanced surfactants, and stable oil-in-water emulsions has shown the ability to significantly reduce viscosity, enabling stable transport even under subzero temperatures (Al-Maamari et al., 2015; Rushd et al., 2023). Simultaneously, advances in drag-reduction strategies—including biomimetic surface treatments, novel drag-reducing agents, and boundary layer control methods—contribute to minimizing hydraulic losses and enhancing throughput (Rushd et al., 2023; Kumar and Mahto, 2017).

Thermal management technologies are another area of progress. Emulsification techniques combined with advanced thermal insulation, such as photothermal coatings, help maintain flow stability and lower heating requirements in cold climates (Jing et al., 2020; Lyu et al., 2019). These solutions not only address technical challenges but also support sustainable pipeline operations by reducing environmental impact and energy demand (Lyu et al., 2019; Al-Maamari et al., 2015).

Overall, recent advances reflect a shift toward integrated strategies that combine rheological modeling, computational methods, chemical treatments, and thermal technologies. Their application in industry demonstrates that modern heavy oil transportation research is increasingly focused on balancing operational efficiency, cost-effectiveness, and environmental responsibility.

Despite the wide range of viscosity and rheological models proposed for crude oil mixtures, several fundamental limitations remain unresolved when these models are applied to non-Newtonian heavy oils under cold climate conditions. Most existing correlations account for either temperature effects or compositional variation, but rarely address their coupled influence on non-Newtonian flow behavior. In particular, classical Arrhenius- and Refutas-type models neglect shear-dependent effects altogether, while conventional power-law and Herschel–Bulkley formulations typically assume constant rheological parameters, which is physically inconsistent for systems undergoing thermal and compositional changes.

Furthermore, commonly used mixing rules implicitly assume ideal or near-ideal blending behavior and therefore fail to capture nonlinear inter-component interactions, especially in mixtures exhibiting abrupt rheological transitions associated with asphaltene aggregation and paraffin crystallization at low temperatures. As a result, these models often lose predictive accuracy under cold climate conditions and at high heavy-oil concentrations, precisely where reliable rheological prediction is most critical for flow assurance.

The present work addresses these limitations by proposing a physically based rheological model in which both the consistency coefficient and the flow behavior index are expressed as explicit functions of temperature and mixture composition, augmented by interaction parameters that account for nonlinear structural effects. Unlike conventional approaches, the proposed formulation enables continuous prediction of rheological behavior across wide temperature and concentration ranges while preserving correct limiting behavior for pure components. This provides a unified framework for predicting non-Newtonian flow behavior of crude oil mixtures under cold climate conditions.

The main contributions of this study are threefold:

Development of a rheological model that simultaneously accounts for temperature, composition, and non-Newtonian effects in crude oil mixtures.

The proposed rheological model was developed following a stepwise methodology that integrates physical reasoning, experimental observations, and empirical parameter identification. As a starting point, the Ostwald–de Waele (power-law) model was selected as the base formulation due to its widespread use in describing shear-thinning behavior of non-Newtonian crude oils. However, in contrast to classical implementations with constant rheological parameters, the present approach assumes that both the consistency coefficient and the flow behavior index vary systematically with temperature and mixture composition.

In the first stage, temperature effects were incorporated by introducing exponential temperature-dependent expressions for the consistency coefficient and the flow behavior index. This reflects the physical weakening of intermolecular interactions, including asphaltene and paraffin associations, with increasing temperature. In the second stage, concentration effects were introduced by expressing these rheological parameters as functions of the heavy crude oil fraction in the mixture. To account for non-ideal blending behavior, inter-component interaction terms were included, enabling the model to capture nonlinear and asymmetric effects observed during dilution.

In the final stage, additional correction parameters were introduced to describe experimentally observed abrupt rheological transitions occurring at low temperatures and high heavy-oil concentrations. These parameters represent concentration-induced structural instabilities and allow the model to reproduce sharp changes in viscosity and shear stress that cannot be captured by smooth interpolation alone. The resulting formulation provides a unified expression for shear stress and apparent viscosity as simultaneous functions of temperature, composition, and shear rate.

The experimental component of the methodology was designed to provide comprehensive rheological data for model calibration and validation. Binary mixtures of high-viscosity Severo-Komsomolskoye crude oil and low-viscosity Vankor crude oil were prepared over a wide range of volumetric ratios, covering the full composition interval from pure light to pure heavy crude. Rheological measurements were performed using a rotational rheometer under controlled thermal conditions.

Experiments were conducted over a temperature range of 5 °C–60 °C and shear rates from 1 to 300 s^-1, corresponding to operating conditions relevant to pipeline transport in cold regions. A full factorial multi-level experimental design was employed to systematically investigate the combined influence of temperature and mixture composition. Detailed descriptions of materials, sample preparation, measurement procedures, and experimental design are provided in the Materials and Methods section.

Model parameters were identified through systematic processing of the experimental data. For each temperature and mixture composition, shear stress–shear rate data were analyzed in logarithmic coordinates to determine the consistency coefficient and flow behavior index from linear regression. These experimentally derived parameters formed the basis for subsequent calibration of the temperature- and concentration-dependent model functions.

Empirical coefficients governing temperature dependence, concentration effects, and inter-component interactions were determined by minimizing the deviation between experimentally measured and model-predicted shear stress values over the entire dataset. To evaluate predictive robustness and avoid overfitting, the experimental data were divided into calibration and validation subsets, with approximately 70% of the data used for parameter estimation and 30% reserved for independent validation.

Model performance was assessed using statistical error metrics, including the mean relative deviation and root mean square error. The resulting calibrated model was then used to predict shear stress and viscosity across the full range of investigated temperatures, compositions, and shear rates, enabling direct comparison with experimental observations in the Results and Discussion section.

The proposed rheological model is applicable within the range of conditions covered by the experimental program, namely temperatures from 5 °C to 60 °C, heavy crude oil concentrations from 0% to 100%, and shear rates between 1 and 300 s^-1. The model is formulated for steady-state shear flow and was calibrated using laboratory-scale rheological measurements. Outside this range, particularly at temperatures below 5 °C or under transient and viscoelastic flow conditions, additional structural effects may arise that are not explicitly accounted for in the present formulation. Therefore, extrapolation beyond the validated domain should be performed with caution.

The experimental investigation was carried out using crude high-viscosity oil samples obtained from the Severo-Komsomolskoye oil field and low-viscosity samples of the Vankor oil field (Figure 6). To ensure a comprehensive investigation of rheological behavior, binary mixtures of the Severo-Komsomolskoye (high-viscosity) and Vankor (low-viscosity) crudes were prepared in volumetric ratios of 1/9, 2/8, 3/7, 4/6, 5/5, 6/4, 7/3, 8/2, and 9/1. These proportions were selected to provide detailed coverage of the entire compositional range, allowing systematic evaluation of viscosity and shear-dependent flow properties across different blending ratios.

Rheological measurements were performed using a rotational rheometer Rheotest RN 4.1 (Figure 7), which allows precise determination of shear stress, shear rate, and viscosity under controlled thermal conditions. The instrument operates by measuring the torque required to rotate one cylindrical or conical element relative to another, providing direct evaluation of shear-dependent viscosity.

Experiments were conducted in a shear rate range of 1–300 s^-1 and a temperature range of 5 °C–60 °C, corresponding to expected operational conditions in northern pipeline systems. The measurement protocol included three stages (Figure 8):

CR sweep (0 → 300 s^-1), where the shear rate was gradually increased;

Prior to testing, samples were degassed, thermostabilized, and equilibrated for 15–20 min to ensure uniformity. Each test was repeated at least three times to ensure reproducibility.

A full factorial multi-level design was employed to systematically study the influence of two independent factors:

Temperature (X₁): 5 °C–60 °C;

Concentration of Severo-Komsomolskoye crude in mixture (X₂): 0%–100%.

The experimental design was structured according to standard multifactorial planning methods (Huan et al., 2024; Jackson, 1994), including both two-level and three-level factorial designs, which were extended to a full multi-level matrix to improve accuracy. To this end, temperature and high-viscosity oil concentration factors were coded as shown in Table 2.

In total, 88 experimental runs were performed to cover the specified ranges of temperature and concentration (Table 3).

Shear stress (τ) and dynamic viscosity (η) were determined for each experimental point using the relationship as Equation 21 η = τ γ ˙ (21) where γ ˙ is the applied shear rate.

Data were processed to eliminate measurement artifacts (e.g., local heating effects at high shear rates). Averaged results were used to construct flow curves τ γ ˙ and viscosity curves η γ ˙ at different temperatures and mixture compositions. The results were further analyzed to identify Newtonian and non-Newtonian regimes, yield stress behavior, and temperature–concentration effects.

The empirical parameters α and β were calibrated simultaneously with the remaining model coefficients using the full experimental dataset (88 runs). Their values were selected to ensure accurate reproduction of experimentally observed inflection points and regions of rapid rheological transition, which could not be captured by conventional power-law formulations.

At the first stage of data processing, only the experimental points that were statistically significant and representative of steady-state flow conditions were selected for further analysis. All anomalous points associated with measurement fluctuations or instrument stabilization were excluded.

For each temperature and concentration, plots of the dependence ln τ = ln K + n · ⁡ ln γ ˙ were constructed. From the slope and intercept of these linear plots, the flow behavior index n and the consistency coefficient K were determined for each experimental condition.

Subsequently, the obtained K and n values were organized into tables showing their dependence on both temperature and concentration. Separate datasets were compiled for K T , K C , n T , and n C , which served as the basis for parameter estimation of the proposed model.

Initial approximate values of all empirical coefficients in the governing equations were then assigned. By iteratively adjusting these empirical parameters, the deviation between the experimentally determined and model-calculated values of the consistency coefficient K and the flow index n was minimized. The optimal set of parameters was obtained when the total mean-square deviation between predicted and experimental data reached its minimum, ensuring the best fit of the model to the experimental results.

The empirical parameters of the proposed rheological model were determined using an iterative nonlinear regression procedure. As an objective function, the mean squared deviation between experimentally measured and model-predicted shear stress values was minimized over the calibration dataset as Equation 22 J = 1 N ∑ i = 1 N τ i exp − τ i c a l c 2 (22) where N is the number of experimental data points. Initial parameter estimates were obtained from linear regressions of the logarithmic shear stress–shear rate data and from physically reasonable bounds inferred from the experimental trends. Model parameters were then updated iteratively to reduce the objective function value.

Convergence was assumed when the relative change in the objective function between successive iterations fell below 10^–6, and when further parameter updates produced changes in the predicted shear stress smaller than the experimental uncertainty. The optimization process was additionally constrained to ensure physically meaningful parameter values, including positivity of the consistency coefficient and monotonic trends of rheological parameters with temperature. The final parameter set corresponded to the global minimum of the objective function within the explored parameter space.

To assess the predictive robustness of the proposed model, the experimental dataset was partitioned into two subsets. Approximately 70% of the data were used for parameter calibration, while the remaining 30% served as an independent validation set. The partitioning was performed to ensure representative coverage of the full temperature and concentration ranges in both subsets.

The predictive accuracy obtained for the validation subset was comparable to that of the calibration data, indicating stable model behavior and the absence of overfitting.

After determining the main parameters K and n , the shear stress was calculated at a constant shear rate using the proposed rheological model. The obtained values were plotted together with the experimental data for comparison. The resulting graph, presented in Figure 9, demonstrates a close agreement between the experimental and calculated shear stress values, confirming the validity of the proposed model under the studied conditions.

The results obtained for the shear stress at different temperatures and concentrations of heavy crude oil are summarized in Table 4. The third column presents experimentally measured values, while the fourth column corresponds to those calculated using the proposed rheological model.

A quantitative comparison between the experimental and model-predicted shear stress values demonstrates a high level of agreement. The relative deviation δ between experimental τ exp and calculated τ calc values was evaluated using the expression as Equation 23 δ = τ exp − τ calc τ exp · 100 % (23)

The mean relative deviation across all data points did not exceed 8.7%, while the root mean square error (RMSE) was less than 0.95 Pa, which indicates excellent predictive accuracy for practical engineering applications.

In addition to global performance metrics, the distribution of relative deviations was analyzed. The residuals exhibited an approximately symmetric distribution centered near zero, indicating the absence of systematic bias across the investigated operating conditions.

Assuming a normal distribution of residuals, a 95% confidence interval for the predicted shear stress was estimated based on the RMSE. For most temperatures and mixture compositions, the confidence bounds remained within ±1.9 Pa, confirming the statistical stability of the model predictions.

A comparison of the data demonstrates that the model reproduces the experimental results with high accuracy across the entire range of studied temperatures (5 °C–60 °C) and concentrations (0–1). The deviations between the experimental and calculated shear stress values do not exceed 10%–12% for most conditions, confirming the reliability of the adopted parameterization and the physical adequacy of the model.

At higher temperatures (50 °C–60 °C), the experimental and calculated curves nearly coincide, indicating that the model correctly captures the transition of the mixture toward Newtonian flow. In this temperature region, the weakening of molecular interactions leads to reduced structural viscosity, and the flow behavior index n approaches unity.

A systematic analysis of prediction errors indicates that the largest deviations occur at temperatures below 15 °C and heavy crude concentrations exceeding 0,7. These conditions represent the boundary of the validated applicability range, where strong structural heterogeneities and time-dependent effects become increasingly significant. In this regime, the relative deviation increases due to the enhanced contribution of temperature-sensitive structural phenomena, such as asphaltene aggregation and paraffin crystallization.

The increased deviations observed at temperatures below approximately 15 °C can be attributed to physicochemical phenomena that are only partially captured by the steady-state power-law framework. At low temperatures, heavy crude oil mixtures exhibit enhanced structural heterogeneity due to intensified asphaltene aggregation, paraffin crystallization, and the formation of transient network structures. These processes introduce time-dependent and history-dependent rheological effects, including thixotropy and weak viscoelasticity, which are not explicitly represented in the present formulation.

In addition, low-temperature conditions amplify the sensitivity of rheological properties to minor compositional fluctuations and experimental uncertainties, particularly at high heavy-oil concentrations. As a result, localized deviations between experimental and predicted shear stress values become more pronounced. Future extensions of the model may incorporate structural kinetics, viscoelastic constitutive terms, or temperature-dependent yield stress formulations to better capture these effects. Nevertheless, within the validated temperature and concentration range, the current model provides stable and physically consistent predictions suitable for engineering calculations.

The dependence of shear stress on concentration at a constant shear rate (Figure 4) shows a clearly nonlinear character. At low concentrations ( C < 0.4 ), the addition of a light component causes a rapid decrease in shear stress—a typical dilution effect. However, at higher concentrations ( C > 0.7 ), the effect becomes weaker, and the mixture behaves similarly to a structured non-Newtonian fluid. The experimental and model-predicted surfaces coincide closely, confirming that the proposed Equations 8–10 adequately describe the joint influence of thermal and compositional factors on the rheological properties of binary crude oil mixtures.

Overall, the results validate the applicability of the developed model for predicting shear stress and viscosity under a wide range of operating conditions. The high level of agreement between experimental and calculated data demonstrates that the inclusion of inter-component interaction coefficients significantly improves the accuracy of rheological predictions compared with classical models.

A qualitative sensitivity assessment indicates that model predictions at low temperatures are primarily governed by the consistency coefficient parameters, whereas at higher temperatures and lower heavy oil concentrations, the flow behavior index becomes the dominant factor.

To demonstrate the practical use of the proposed rheological model, its predictions for the same experimental conditions were directly compared with those obtained using widely adopted Arrhenius, Refutas, and classical power-law models.

To quantitatively assess the performance of the proposed rheological model, its predictive accuracy was compared with several widely used viscosity and rheological correlations reported in the literature. Specifically, the Arrhenius temperature–viscosity model, the Refutas viscosity blending rule, and the classical power-law model with constant parameters were applied to the same experimental dataset.

All reference models were calibrated using the same experimental data and evaluated under identical temperature and concentration conditions to ensure a fair comparison (Table 5). The comparison is performed using the same temperature, composition, and shear rate conditions as those employed in the validation example presented in Figure 9 and Table 4.

The comparison demonstrates that conventional viscosity correlations exhibit substantially higher prediction errors when applied to non-Newtonian crude oil mixtures under varying temperature and compositional conditions. Models accounting for only one governing factor—temperature, concentration, or shear dependence—fail to reproduce the coupled effects observed experimentally.

In contrast, the proposed model explicitly incorporates temperature-dependent non-Newtonian behavior and inter-component interactions, resulting in a significantly lower prediction error across the entire operating domain. This confirms the practical advantage of the model for pipeline transport applications in cold climate conditions.

The superior predictive performance of the proposed model is consistent with the qualitative sensitivity trends discussed above, further confirming its robustness across the investigated operating conditions.