Assessing Rheology Effects and Pore Space Complexity in Polymer Flow Through Porous Media: A Pore‐Scale Simulation Study

Non‐Newtonian fluid flow within porous media, exemplified by polymer remediation of contaminated groundwater/aquifer systems, presents complex challenges due to the fluids' complex rheological behavior within 3D tortuous pore structures. This paper introduces a pore‐scale flow simulator based on the OpenFOAM open‐source library, designed to model shear‐thinning flow within porous media. Leveraging this developed solver, extensive pore‐scale flow simulations were conducted on μ‐CT images of various real and synthetic porous media with varying complexity for both power‐law and Cross‐fluid models. We focused on the macroscale‐averaged deviation between bulk viscosity and the in‐situ viscosity, commonly denoted by a shift factor. We provided an in‐depth evaluation of the shift factor's dependency on the fluid's rheological attributes and the rock's pore space complexity. The least‐squares fitted values of the shift factor fell in the range of 1.6–9.5. Notably, the most pronounced shift factor emerged for extreme flow behavior indices. Our findings highlight not just the critical role of rheological parameters, but also demonstrate how the shift factor fluctuates based on tortuosity, characteristic pore length, and the cementation exponent. In particular, less porous/permeable systems with smaller characteristic pore lengths exhibited larger shift factors due to higher variations of shear rate and local viscosity in narrower flow paths. Additionally, the shift factor increased as rock became more tortuous and heterogeneous. The introduced pore‐scale simulation proves instrumental in determining the macroscopic averaged shift factor. This, in consequence, is vital for precise computations of viscosity and pressure drop when dealing with non‐Newtonian fluid flow in porous media.


Introduction
The flow of polymeric liquids within porous materials has significant relevance across various applications, including polymer flooding in hydrocarbon reservoirs (Kamal et al., 2015;Mohsenatabar Firozjaii & Saghafi, 2020), remediation of contaminated aquifer systems (Gastone et al., 2014;Skauge et al., 2018), treatment of contaminated surface site formations (Zhong et al., 2011), enhancing the hydraulic fracturing performance (Hauswirth & Miller, 2014), and fiber molding (Li et al., 2016).Polymeric solutions stand out as a notable subset of non-Newtonian fluids.Understanding the interplay between pressure drop and flow rate for diverse fluid compositions and velocity spectra is vital, especially when predicting groundwater and crude oil extraction rates from subsurface reservoirs or designing industrial operations that require the flow and transportation of non-Newtonian fluids through porous matrices.However, due to complex interactions, the flow dynamics of polymer solutions in porous media remain somewhat elusive, largely attributed to the intricate design and morphology of pore spaces coupled with the fluid behavior.These factors raise difficulty in describing flow physics, consequently hindering accurate predictions of macroscopic pressure drops.
Contrasting Newtonian fluids, the viscosity of non-Newtonian fluids shifts based on the exerted force -more precisely, the viscosity solely depends on the shear rate during steady, simple shear flow.Thus, stress and strain rate exhibit a nonlinear relationship (Bird et al., 1987).Broadly, non-Newtonian fluids can be segmented into two categories (Muljadi et al., 2016); the first encompasses rheopectic and thixotropic fluids where viscosity is timedependent concerning the applied force, and the second consists of dilatant (or shear-thickening), pseudoplastic (or shear-thinning), and Bingham fluids where viscosity remains time-invariant.Shear-thinning fluids, characterized by a decreasing viscosity with escalating shear rates, are perhaps the most ubiquitous non-Newtonian fluids, finding extensive practical applications.The challenge of modeling non-Newtonian flow is amplified when confined to the intricate pore spaces of porous materials.For cases where the correlation between fluid velocity and the hydraulic gradient is linear, as seen with Newtonian fluid dynamics, Darcy's law is the applied phenomenological equation to describe steady-state flow through porous media under low Reynolds number (Re) conditions.Yet, as the Re increases, nonlinearities emerge.The roots of such nonlinear patterns can be traced back to factors like microscale drag forces (Zhang et al., 2019), inertial effects (Hassanizadeh & Gray, 1987), and low Re turbulence (Levy et al., 1999).A common approach adopted in the literature for the laminar flow of non-Newtonian fluids through a porous medium is to use the Darcy equation with all non-Newtonian effects combining into a single viscosity parameter termed the porous media viscosity (Dybbs & Edwards, 1984;Sorbie, 1991).The porous medium viscosity is a function of the so-called porous medium shear rate, which is typically defined based on the superficial fluid velocity and the characteristics of the medium (Eberhard et al., 2019;Sorbie, 1991).
Experimental investigations have highlighted that the behavior of the viscosity of the bulk non-Newtonian fluid differs from that within porous media (Tosco et al., 2013).In order to account for such discrepancies, researchers have introduced a shift or scaling factor (Hauswirth et al., 2020;Sorbie et al., 1989;Willhite & Uhl, 1988).This factor has proven instrumental in fitting experimental data (Sheng, 2010).Previous research has shown that the shift factor is influenced by both the rheological properties of the fluid and the intrinsic properties of the porous medium, such as tortuosity (Tosco et al., 2013;Zhang et al., 2019).Values for this shift factor, as documented by several researchers, typically range between 1 and 15 (Hauswirth et al., 2020;Jithin et al., 2018;Tosco et al., 2013).On macroscopic scales, several studies have proposed semiempirical equations for the shift factor when considering non-Newtonian fluid models-predominantly the power-law model-as they move through porous materials, based on the capillary-bundle models of porous media (Christopher & Middleman, 1965;Eberhard et al., 2019;Hirasaki & Pope, 1974;Sorbie et al., 1989;Teeuw & Hesselink, 1980).Knowing the bulk rheological characteristics of the fluid, the permeability of the porous media, and the shift factor is required for the macroscale description of the dependence of pressure drops on flow velocity.These semiempirical models, however, do not provide insights into the relationship between the shift factor (as a macroscale-averaged parameter) and the microscale parameters such as pore geometry and topology.
In a pioneering work, Sorbie et al. (1989) utilized network modeling with a 2D capillary array and the Carreau fluid model to simulate polymer flow in porous media.They obtained shift factors ranging from 1.2 to 1.67 for the network and pointed out larger ones for more heterogeneous networks.Tosco et al. (2013) used SEM images of natural unpacked sand grains (porosity range of 0.4-0.5) to perform numerical flow simulations for both Newtonian and non-Newtonian fluids using CFD code ANSYS Fluent based on a finite volume discretization.By testing various shear-thinning fluid models (Cross, Ellis, and Carreau models), they evaluated the interplay of porous media and fluid properties on equivalent viscosity and macroscopic pressure gradient and found that the shift factor and the inertial coefficient were dependent on both the porous medium and fluid properties.However, their study was primarily for 2D image slices, putting aside the 3D pore-space connectivity.Jithin et al. (2018) employed the lattice-Boltzmann method to simulate the Carreau-Yasuda fluid flow within both idealized and stochastically reconstructed pore structures.They assessed the impacts of flow and geometric parameters and found that the shift factor was influenced by multiple parameters, including porosity, permeability, tortuosity, and percolation threshold.Zhang et al. (2019) conducted a 3D microscale flow simulation using OpenFOAM in a rough fracture for Cross shear-thinning fluids with varying rheological parameters and over a range of flow regimes.They showed that the shift factor was intertwined with fluid rheology and fracture aperture variation.
Further, Rodríguez de Castro and Agnaou (2019) examined the Herschel-Bulkley and the Carreau fluid flow in various media, highlighting the interaction of yield stress and pore morphology on the shift factor's relationship with flow rate.In particular, Rodríguez de Castro and Radilla (2017b) conducted experiments with xanthan gum solutions as shear-thinning fluids to understand how fluid rheology impacts flow rate and pressure drop in glass sphere packs.In a more recent exploration, Hauswirth et al. (2020) used the LBM approach to simulate Crossfluid flow through synthetically created porous media.They found good agreement between the simulation and experimental results and estimated the values of the shift factor for different sphere-pack models.This study also contrasted the performance of the LBM solver against the OpenFOAM simulator.Yet, amidst these strides, a discernible gap exists: the dynamics of complex non-Newtonian fluids flowing through real porous structures, especially sedimentary rocks, remain largely underexplored.Our current study seeks to bridge this gap by examining complex rheological models, like the Cross-fluid model, along with the simple power-law model, within real-world sandstone and carbonate samples.
Modeling non-Newtonian flow in porous media presents inherent complexities where not only the individual properties of the fluid and the porous medium are essential, but there is also the need to describe the shift factor-a shared property of both elements.Historically, correlations for the shift factor have relied heavily on the capillary bundle model of porous media.However, its application to complex real-world porous rocks remains a matter of debate, underscoring the need for more exhaustive studies.Furthermore, most semi-empirical models, rooted in the assumption of a constant tortuosity, deduce that the shift factor's variation arises solely from the fluid's rheological properties, like the flow behavior index.This is further emphasized by the predominant literature focus on the power-law rheological model, which assumes a straightforward linear correlation between non-Newtonian viscosity and shear rate in a log-log scale.Such an assumption may be inaccurate for many non-Newtonian fluids, such as polymeric solutions, which often deviate from the power-law behavior.Consequently, there is a need for more accurate and realistic models.While numerous research studies over the past decade have aimed to determine the shift factor via microscale flow simulations, the predominant focus has been on idealized 2D/3D structures.This leaves a knowledge gap regarding the applicability of these findings to real 3D porous materials, highlighting the need to critically assess the shift factor in real porous microstructures.The work presented in this study should be understood concerning the lack of sufficient and accurate research on the shift factor for real 3D microstructures, and it attempts to provide insight into the interplay between the shift factor, pore structures, and fluid characteristics.To achieve this, we have performed a large number of 3D porescale flow simulations, employing both Newtonian and non-Newtonian fluids, to probe the shift factor's dependence on fluid rheology and rock pore space characteristics.A numerical solver based on the OpenFOAM finite volume library has been developed to map the flow of non-Newtonian fluids-specifically the power-law and Cross models-in porous media and compute the shift factor.To validate our approach, we have leveraged μ-CT images of various consolidated and unconsolidated porous media.The study also examines the implications of pore space complexities on the shift factor, considering attributes such as tortuosity, characteristic pore length, and cementation exponent.Our simulations span a wide shear rate spectrum in the low (creeping) flow regime.

Theoretical Background
In this section, we first introduce the local governing and constitutive equations for non-Newtonian fluid flow through the pore space of porous materials.Then, we present the numerical techniques used for solving the governing equations and calculating the required parameters from pore-scale flow simulations.

Rheological Behavior of Polymeric Fluids
The flow resistance to the shear rate of simple fluids is characterized by a physical property known as viscosity, which is independent of the shear rate and is only a function of temperature and pressure (Bird et al., 2002).Unlike

10.1029/2023WR036125
Newtonian fluids, a constant viscosity parameter cannot describe the rheological behavior of complex fluids, such as polymeric liquids.A simple common approach to describe the non-Newtonian viscosity, neglecting timedependent effects, is to use the generalized Newtonian fluid model (Bird et al., 2002): where τ is the shear stress, that is, the shearing force per unit area, γ is the shear rate, that is, the rate of deformation or velocity gradient, and µ( γ) is the non-Newtonian viscosity, which is a function of shear rate.Several empirical non-Newtonian viscosity functions, also known as the constitutive rheological equations, have been proposed, of which the simplest one is the two-parameter power-law model (Bird et al., 1987): where K is the flow consistency index, and n is the power-law index (or the flow behavior index).The power-law expression describes the linear portion of the log-log plot of the viscosity versus shear rate curves, known as the viscosity (or flow) curves.In this study, we focus on the cases where n is less than one, that is, the cases where the viscosity decreases with increasing shear rate.This rheological behavior, known as shear-thinning or pseudoplastic, is the most common behavior on non-Newtonian fluids, such as a majority of polymer solutions (Darby et al., 2017).
In practice, however, more accurate constitutive rheology models are required for polymeric fluids (Darby et al., 2017).Hauswirth et al. (2020) reported some well-known examples of non-Newtonian rheological models.
In this study, we consider the four-parameter Cross-fluid model given by Sochi (2010): where μ 0 is the zero-shear-rate viscosity, μ ∞ is the infinite-shear-rate viscosity, and m is a time constant (Tosco et al., 2013).It is worth noting that the exponent n in the power-law model has the same meaning as the n in the Cross equation.

Governing Equations for Flow of Non-Newtonian Fluids in Porous Media
In the recent two decades, direct modeling has attracted much attention for computing single-and multiphase flow in porous media (Blunt et al., 2013).In this approach, the fundamental equations of fluid mechanics, namely the continuity and momentum equations, are directly applied to the pore space of porous media derived from binarized digitized images.For the flow of non-Newtonian incompressible fluids through the pore space of porous materials at steady-state conditions, the governing equations are the continuity and momentum equations given by Bird et al. (1987): where v is the velocity vector, p is pressure, and ρ is fluid density.The common boundary conditions (BCs) used in the literature are the no-slip condition, where the fluid velocity is set to zero on grain boundaries.To determine velocity distribution from Equations 4 and 5 with proper BCs, one must use proper models for the shear stress in terms of velocity gradients (shear rate) and fluid properties.The generalized Newtonian fluid model (Equation 1) is usually used because of its simplicity.In addition, it is required to use empirical expressions for the non-Newtonian viscosity, for example, Equation 2 or 3.
Once pressure and velocity fields are computed from pore-scale flow simulations, obtained by solving Equations 4 and 5 on the pore space region of a porous medium, one can determine the macroscale-averaged parameters (permeability and shift factor) discussed in the following.Permeability is defined by Darcy's law, which describes a linear relationship between pressure drop and superficial velocity for the laminar flow of Newtonian Water Resources Research 10.1029/2023WR036125 fluid through a fully saturated porous medium (Darcy, 1856).The one-dimensional Darcy's equation along the xdirection is given by: where u 0 is the x-component of superficial velocity, k is the permeability of the porous medium, μ is the fluid viscosity, and ∂p/∂x is the pressure gradient along the x-axis.When the velocity of the fluid increases, the flow regime will eventually deviate from Darcy's law due to the inertial effects.The velocity corresponding to the onset of the non-Darcy flow (nonlinear relationship between the pressure gradient and the superficial velocity) is commonly referred to as the critical velocity.
For the laminar flow of non-Newtonian fluids in porous media, a modified Darcy's law can be adapted in combination with the generalized Newtonian fluid model (Equation 1).The adopted approach assumes that the same equations used for Newtonian flow can be utilized for non-Newtonian flow, provided that the constant viscosity μ is replaced by μ pm , the porous medium viscosity, which depends on the properties of the porous medium and the flowing fluid (Bird et al., 1987;Sorbie, 1991).For the case of non-Newtonian flow, Equation 6is written as: where γpm is the porous medium shear rate, defined as the average shear rate that the fluid experiences when flowing through the pore space of the porous medium (Pearson & Tardy, 2002).The porous medium viscosity is used to describe the observed macroscopic rheology of the non-Newtonian fluid in a porous medium (Tosco et al., 2013).The rheological constitutive models are modified to incorporate the effects of the porous medium on flow.For example, the power-law (Equation 2) and Cross fluid (Equation 3) models for a non-Newtonian fluid flow in a porous medium are modified as:

The Capillary-Bundle Model Approach
The bundle-of-capillary-tubes model is one of the most common approaches widely used to model flow and transport in porous media.For the flow of Newtonian fluid in a bundle of capillary tubes with different radii, the equivalent or hydraulic radius, R eq , is defined from a combination of Darcy's and Poiseuille's laws as follows (Eberhard et al., 2019;Rodríguez de Castro & Agnaou, 2019): where ϕ is porosity and τ is tortuosity, defined as the ratio of the actual flow path inside a porous medium to the length of the porous medium (Ghanbarian et al., 2013).Theoretically, for Newtonian fluid flow through a straight bundle of capillaries of uniform diameter, the shear rate at the wall of the capillary tube (denoted by γw,NF ) is given by Chauveteau and Zaitoun (1981), Sorbie et al. (1989), and Willhite and Uhl (1988): For the flow of non-Newtonian fluids through a porous medium, the common approach in the literature is to define the porous medium shear in terms of fluid superficial velocity and properties of the medium (Lopez et al., 2003;Perrin et al., 2006;Sorbie, 1991;Sorbie & Huang, 1991;Valvatne et al., 2005), given by: Water Resources Research where α is a constant empirical value known as a shift or scaling factor.Fundamentally, the shift factor is related to the deviations observed between the viscosity curves (rheograms) of the bulk fluid and those of the fluid within the porous medium (Zhang et al., 2019).The shift factor is basically an averaged macroscopic property, which, unlike permeability, is a function of both the rheological properties of the fluid and the intrinsic properties of the porous medium, for example, the tortuosity (Lopez et al., 2003;Pearson & Tardy, 2002;Sorbie et al., 1989;Tosco et al., 2013).
Several authors have developed analytical and semi-empirical equations for shear rate for the simplest case of laminar flow of power-law fluids through simplified geometries such as tubes, slits, bundles of capillary tubes, and more complicated porous media such as packed beds (Berg & van Wunnik, 2017;Sorbie et al., 1989;Tosco et al., 2013;Zhang et al., 2019).For example, the corresponding equation to Equation 12for power law flow is given by Bird et al. (2002) and Willhite and Uhl (1988): where γw,PL is the shear rate at the wall of the straight capillary tube for the power-law flow.Based on the capillary bundle model of porous media and the Blake-Kozeny equation, the following general model can be proposed for power-law fluid flow in a porous medium: where C and D are constants whose values were reported differently in the literature, as shown in Table 1.The general expression of shift factor α derived based on the capillary bundle model is given by: Table 1 shows the expressions for α based on different investigators.Based on this model, the shift factor is only a function of the power-law index, as evident from Table 1.The derivation of expressions given in Table 1 was based on some preliminary assumptions.First, an equivalent capillary tube radius (Equation 10) was used.Second, to account for the tortuous flow paths, a constant value was replaced for the tortuosity.For example, the tortuosity of the flow paths through the packed beds of spherical particles has been widely reported to be 25/12 (Christopher & Middleman, 1965;Willhite & Uhl, 1988).

Characterization of Pore Structure Complexity
In this study, we consider the following characteristics of a porous medium as indicators of the medium's pore complexity: Hirasaki and Pope (1974 Teeuw and Hesselink (1980) Water Resources Research 10.1029/2023WR036125 The square root of permeability (√k) is a common characteristic length scale of a porous medium (Kececioglu & Jiang, 1994;Ward, 1964).It can be argued that the more microporous the porous medium (usually carbonate rocks), the less the permeability and hence the smaller characteristic length.The product of the formation factor and porosity, F .ϕ, is a measure of the electrical tortuosity of a porous medium (Clennell, 1997;Ghanbarian et al., 2013;Torquato, 2002).A physical meaning can be associated with F .ϕ so that for a bundle of straight capillary tube model, F .ϕ = 1.The more tortuous the pore space, the greater the value of F .ϕ.The cementation exponent, m c , is determined from the Archi equation (Clennell, 1997): For a bundle of straight capillary tubes, m c = 1 (Glover, 2009).For real rocks, the cementation exponent varies for different rock types and practically ranges between 1.5 and 4 (Archie, 1942).For a simple case of spherical rock grains, Sen et al. (1981) found from the analytical investigation the cementation exponent is equal to 1.5.For consolidated rocks, m c is usually around 2 (Torquato, 2002).In general, more heterogeneous porous rocks have higher values of m c .Glover (2009) argued that the pore space connectivity depends on the porosity and the cementation exponent.

Conditions of Laminar Non-Newtonian Flow in Porous Media
This study considers the laminar flow of non-Newtonian fluids through porous media.A common approach in the literature is to use Reynolds numbers to identify the creeping (laminar) flow and mark the beginning of the inertial flow regime (Koch & Hill, 2001;Sederman et al., 1998).For simple grain-based models of porous materials, the onset of inertial effects for the flow of Newtonian fluids is usually expressed in terms of Reynolds numbers based on the mean grain diameter as a characteristic length (Arbabi & Sahimi, 2024).For more complex systems, such as natural rock materials, a common approach is to define the Reynolds number as follows (Beavers & Sparrow, 1969;Muljadi et al., 2016;Wood et al., 2020) In Equation 17, μ is the fluid viscosity and the square root of permeability is used as the characteristic pore length.
The problem is much more complicated when dealing with non-Newtonian flow through porous materials.
Several approaches have been used to define the critical Reynolds numbers to indicate the end of the creeping flow (Chhabra, 2006) in non-Newtonian flow.Over the past several decades, a considerable amount of experimental and theoretical research has focused on determining the pressure drop of non-Newtonian fluid flow through porous media and estimating the limits of creeping flow regime mainly by assuming the power-law model (Comiti et al., 2000;Morais et al., 2009;Sabiri & Comiti, 1995).In this study, we adopted a unified definition for the Reynolds number, denoted by Re Non-Newtonian , which is applicable to both the power-law and Cross-fluid models in identifying the onset of inertial flow effects.In this approach, Re Non-Newtonian is obtained by modifying Equation 17to incorporate the porous medium viscosity instead of Newtonian viscosity, as given by Equation 18: It is noted that similar approaches were used in the literature to define the Reynolds number for non-Newtonian flow in porous media (Shende et al., 2022;Zhang et al., 2019).

Description of μ-CT Images and Pore-Space Complexity
The  2 reports the image dimension and voxel size used in flow simulations.The single-phase Newtonian fluid flow was first simulated to obtain permeability (k) and electrical formation factor (F), as given in Table 2.
Only the values of k and F along the x-axis were considered in this paper.As evident from Figure 1, the samples display diverse levels of complexity in their pore space.Table 2 reports the values of the indicators used in this study to quantify the complexity of the samples' pore space.It is noted that the values of the cementation exponent were calculated using the Archie equation.

Pore-Scale Flow Simulations and Solver Development
We used the poreFoam software package (https://github.com/ImperialCollegeLondon/poreFoam-singlePhase)developed by Raeini et al. (2012) to directly solve Equations 4 and 5 with proper BCs on the voxelized pore space extracted from μ-CT images of porous media.The poreFoam numerical solver is an open-source CFD software based on the OpenFOAM finite volume library (Shams et al., 2018).The continuity and momentum equations are coupled using the PIMPLE algorithm-a combination of PISO and SIMPLE algorithms.The time differential, ∂(…)/∂t, is solved with an implicit, first-order accurate Euler scheme.The first-order accuracy in time is deemed sufficient given that we are interested only in the steady-state solution (Angrand et al., 1985).The divergence operator, ∇.(…), is discretized with a Gauss scheme and interpolated using a second-order accurate self-filtered central difference scheme (Hafez & Kwak, 2003;Jasak, 1996).
This study modified the poreFoam single-phase solver to model the non-Newtonian flow in porous media.In all simulations, the BCs at the pore-solid interface were set to be a no-slip (zero normal and tangential velocity).In addition, constant pressure BCs at the inlet and outlet faces of the images were used, whereas no-slip BCs were applied on the remaining faces.The criterion used for the steady-state convergence is ε ≤ 10 6 where  In Equation 19, u i and u i+1 are velocities at time steps i and i + 1, respectively.Shear thinning fluids properties for power-law and Cross-fluid models were set in the transport properties.Special attention was paid to the convergence condition using the Courant (or CFL) number (Weller et al., 1998).Following a sensitivity analysis, we set the maximum allowed Courant number equal to 0.4 for the solution of Newtonian fluids to converge.
Regarding the power-law and Cross-fluid models, we set the Courant number limit equal to 0.0001.In both Newtonian and non-Newtonian cases, the time step was chosen to honor the Courant number limit.
The pore-scale flow simulations were conducted for several pressure gradients across the sample.From simulations, we obtained superficial velocity, u, from the input data of macroscopic pressure drop (Δp) across a domain of length L with the assumption of (∂p/∂x) ≈ (Δp/L).We only considered the data on the x-axis.To determine the shift factor, the following procedure was adopted in this study: 1.The permeability of the porous medium, k, is obtained by Darcy's law (Equation 6) using the superficial velocity and pressure drop data provided from the pore-scale flow simulations for Newtonian fluids.2. Flow simulations of non-Newtonian fluids are performed, and μ pm is determined from Equation 7for various pressure drops and flow rates.The simulations are terminated when a deviation from the decreasing trend in the μ pm -u or μ pm -γpm curve (e.g., see Figures 6 and 7) is observed at high flow rates, indicating reaching the critical velocity at which the inertial flow regime begins.3. The shift factor (α) is calculated in terms of simulation data of μ pm (from step 2), flow velocity (u), sample permeability (k) and porosity (ϕ), and the power-law index, n, by the least-square matching using the functional relationship between μ pm and u in Equation 20 or Equation 21.We used the CurveExpert software (https://www.curveexpert.net/)for the purpose of curve fitting.
Power-law model: In this study, the numerical simulation work was implemented on two workstations.A summary of hardware and simulation times for the developed solver is reported in Table 3.Each simulation was carried out in parallel with several processors, depending on the complexity of the sample.The simulation time was varied depending mainly on the applied pressure gradient.With increasing the pressure gradient, the simulation time became larger.For instance, for non-Newtonian fluid flow through the Estaillades sample (image size of 400 3 voxels) at u = 0.195 m/ s, the most challenging case, the run time was about 10 hr.

Solver Verification: Flow of Non-Newtonian Fluids Through a Sphere Pack
To verify the results of pore-scale flow simulations of power-law flow (K = 0.3 Pa.s with varying n) in a bed of uniformly-sized random spherical particles, we extensively performed simulations of power-law flow on the μ-CT image of a beadpack sample named Finney packing of spheres (Finney & Bernal, 1970;Teeuw & Hesselink, 1980).The beadpack image was retrieved from the Digital Rocks Portal archive (http://www.digitalrocksportal.org/projects/47).Figure 1e shows a 2D slice through the 3D segmented image of the sample with the properties reported in  between the analytical cementation exponent value obtained by Sen et al. (1981) and that obtained by the numerical simulation.
The comparison of pore-scale simulations of power-law flow through the beadpack sample and the semianalytical solutions (Equation 14and Table 1) is shown in Figure 2 (Sochi, 2010).Using black solid lines, we also show the model data (Equation 14) based on the estimated (fitted) shift factor value in Figure 2.
The numerical results in Figure 2 were found to be in good agreement with the semi-analytical models (Equation 14 and Table 1).A close look at Figure 2 14) of the power-law fluid flow in a beadpack with our pore-scale numerical simulation results in terms of porous medium viscosity, μ pm , as a function of the intuitive porous medium shear rate, defined by u/√(kϕ).Note that the value of K (the flow consistency index) was taken to be 0.3 Pa.s n .The values of ϕ and k are reported in Table 2. Water for n = 1, for which the power-law model (Equation 8) reduces to the Newtonian fluid model with K = μ.It is noted that all the models given by Equation 14 reduce to μ pm = K for n = 1.From the simulation, the numerical results for n = 1 were found to be 0.29 Pa.s in the creeping flow regime, which is in excellent agreement with the semi-analytical model.In spite of the fact that all the models in Equation 14give the same porous medium viscosity for n = 1, the values of α are different for various models, as indicated in Table 1.As mentioned before, for the flow of Newtonian fluids in a straight bundle of uniform capillaries, the shear rate at the wall of the capillary tube is given by Equation 11.In this case, the equivalent shift factor is equal to (4/√8) ≈ 1.4142, which is the same as the shift factor value obtained by the T-H model for n = 1 (see Table 1).
Figure 3 compares the calculated shift factor (α) values from pore-scale flow simulations and those from the semianalytical models for different values of the power-law index.For the cases where the power-law indices are equal to 0.2 and 0.4, the numerical values of the shift factor are relatively close to those obtained by the C-M model.This is while for the other cases (n = 0.6 and 0.8), an excellent agreement was observed for the numerical values of α and those calculated by the T-H model.

Model Validation: Comparison to Experimental Data
To evaluate the model performance and identify any discrepancies for improvement, we compared the predictions generated by the model with experimental results.Limited experimental data on the shift factors are accessible for comparison.In this study, experimental data collected from two relevant literature cases (Chauveteau, 1982;Rodríguez de Castro & Radilla, 2017b) were used to validate the predictive capabilities of the developed computational model.The experimental setups involved flooding xanthan polymer solutions through packed beds of glass spheres under various conditions.The Carreau model (Equation 22) was used to calculate the shift factor (Bird et al., 1987):  14and Table 1).
where μ 0 , μ ∞ , and n are similar to those in the Cross-fluid model, and λ is a time constant.In Equation 22, γpm is given by Equation 12. Table 4 reports the numerical values of Carreau's parameters used in the simulations for comparison with the experimental data.We utilized Porous Microstructure Generator (PMG) software (Daniel et al., 2023) to generate beadpack numerical models.The grain diameter and porosity of the generated models closely matched or were identical to those of the physical models utilized in the experiments.It is important to highlight that the shift factor was calculated using the same method described in Section 3.3, which involved performing a least-square fit of the simulation data using Equation 22.
Table 5 presents the computed permeability and shift factor values for each numerical model (porous medium and fluid system) alongside the corresponding experimental values of these parameters.Due to variations in permeability and other crucial factors like tortuosity between the physical and numerical beadpack models, there are discrepancies between the experimental and estimated values of the shift factor, although the trends of changes remain consistent across the permeability range.In order to ensure a reasonable comparison between the experimental and simulation results, the shift factor values are normalized by porosity and permeability, taking the form of α/√(kϕ), as illustrated in Figure 4. We observed a strong agreement between the experimental and simulation findings for the experimental data reported by Chauveteau (1982), where permeability ranged between ∼10 mD and ∼140 D (Figure 4a).For the models featuring very high permeability values (Figure 4b) as utilized by Rodríguez de Castro and Radilla (2017b), we noticed a lack of satisfactory agreement by the explanation that the discrepancies may arise due to the departure from Darcy's law under high permeability conditions and the influence of tortuosity on the shift factor.It is noted that the dashed lines depicted in Figure 4 represent the prediction range at a 90% confidence interval, determined by employing the root mean square of the residual between experimental and simulation data as the measure of uncertainty (Longo et al., 2013).Based on the satisfactory match with the experimental results, the developed solver exhibits promising predictive capabilities for modeling the flow of non-Newtonian fluids through porous materials.

Results and Discussion
Table 6 gives the values of parameters of the power-law and Cross-fluid models used in our pore-scale flow simulations.The numerical values of the parameters of the rheological models were chosen based on the data reported in the literature for the Xanthan gum solutions (Sochi, 2010;Soto Caballero et al., 2016;Tosco et al., 2013;Zhong et al., 2013).Flow simulations for various pressure drops and flow rates were performed for different values of the rheological parameters of the shear-thinning fluid (Table 6).It is worth mentioning that the initial Newtonian flow simulations were performed with a constant viscosity of 0.001 Pa.s (i.e., μ 0 = μ ∞ = 0.001 Pa.s in accordance with the data given in Table 6).

The Onset of Inertial Flow Effects
In this study, the simulations were performed until a critical condition marked the beginning of the inertial flow regime.As outlined in Section 2.5 in detail, Equation 18 was used as an appropriate definition for the Reynolds number for the flow of non-Newtonian fluid through porous media.To indicate the transition from the creeping to inertial flow regimes, we define a pressure drop ratio as:  where (Δp/l) apparent is the apparent pressure drop and (Δp/l) creeping is the pressure drop for the creeping flow regime.below 99% of the overall pressure drop despite the sufficiency of 95% for engineering purposes.As evident from Figure 5, regarding the Newtonian fluid model, the critical Reynolds numbers at which the inertial effects begin took place at 6 × 10 3 and 8.5 × 10 5 for Bentheimer and Estaillades, respectively.Our results are in good agreement with the estimation given by Muljadi et al. (2016).
A comparison between the Newtonian and non-Newtonian flow in Figure 5 shows that there are significant fluctuations before the onset of the inertial flow regime for the non-Newtonian flow, particularly for the Cross-fluid model, so it is challenging to choose a clear cutoff (such as 99% or 95% adopted by Muljadi et al. (2016)) to determine the end of the creeping flow regime.Regarding the power-law curves of Figure 5, we used a cutoff value of 95%, whereas a cutoff of 90% was used for the Cross-fluid model due to larger fluctuations.In this respect, the estimated critical Reynolds numbers for all samples for the non-Newtonian fluids with n = 0.62 are shown in Table 7.
From Figure 5, based on the definitions of the Reynolds numbers given by Equations 18 and 23, the critical Reynolds numbers for the non-Newtonian flow are smaller than those of the Newtonian flow by three to six orders of magnitude, representing the importance of rheological effects inside the pore space of porous media.The critical Reynolds number for the Cross-fluid flow in the Bentheimer sandstone is substantially higher than that in the Estaillades carbonate.However, a similar behavior was not observed for the power-law fluid flow in these two rock types, as evident from Figure 5 and Table 7, highlighting complex interactions between the rheological behavior and the pore space complexity; hence, forecasting faces challenges.
Table 7 shows that the critical Reynolds numbers, Re Non-Newtonian , for the Cross-fluid model are one to three orders of magnitude larger than those for the power-law model.Furthermore, from Table 7, it can be seen that the critical Reynolds numbers for Estaillades are considerably smaller than the other samples, while the synthetic silica A1 has the lowest values of the critical Reynolds number.Consequently, the more complicated rock samples in terms of rock type, characteristic pore length, and tortuosity showed a smaller critical Reynolds number for the non-Newtonian flow.

Effects of the Non-Newtonian Flow Behavior Index
Porous medium viscosity is computed using the modified Darcy's law (Equation 7) for different non-Newtonian fluids, porous samples, and flow rates.The influence of the power-law (flow behavior) index (n) on the porous medium viscosity (μ pm ) versus the intuitive porous medium shear rate in the x-direction, defined by ( γpm = u/ √(kϕ)), is illustrated in Figures 6 and 7 for the power-law and Cross-fluid models, respectively.As evident from the figures, a similar trend was observed for all cases, that is, the μ pm decreased with increasing the flow rate before a critical velocity at which flow converted from Darcy to non-Darcy due to the inertial effects.The figures indicate that the value of porous medium viscosity significantly changes with the flow behavior index.A close look at Figures 6 and 7 reveals some important points, as summarized below.
For both the power-law and Cross fluids, the porous medium viscosity decreases as the shear rate increases due to the shear-thinning effects before it begins to increase at high flow rates due to the inertial effects.Regarding the power-law fluid case, Figure 6 reveals that the fluids with lower power-law indices show a greater change in porous-medium viscosity with changing shear rates.This is similar to the behavior of the power-law bulk viscosity (Frankland, 2012).This is while, according to Figure 7, we see a reverse behavior for the Cross-fluid model where porous-medium viscosity varies significantly with the increase of the power-law index.It is worth mentioning again that for the power-law flow with low power-law indices (n = 0.1 and 0.31), we see a deviation from the linear relationship (laminar region) at extremely low flow rates, probably due to the importance of the interaction between the fluid and the pore walls (Sochi, 2010).Figure 6 shows that the porous medium viscosity decreases by increasing the power-law index up to a certain intuitive shear rate (corresponding to a certain superficial velocity).A reverse trend is observed after that certain shear rate.The behavior of the Cross-fluid model is in the opposite direction compared to the power-law fluid flow through the samples (Figure 7).In other words, the porous medium viscosity increases by increasing the power-law index up to a certain shear rate before further reverse behavior.a Based on the 95% cutoff limit.b Based on the 90% cutoff limit.
Water Resources Research

10.1029/2023WR036125
As outlined in Section 3.3, the shift factor (α) was calculated by matching the functional relationship between μ pm and u, Equations 18 and 19, and the simulation data (i.e., the data given in Figures 6 and 7).Tables 8 and 9 report the calculated values of the shift factor for different values of the power-law index, n, for both the power-law and the Cross-fluid models, respectively.The shift factor data was calculated based on the best fit with the model with a coefficient of determination (R squared) greater than 0.98.For each sample, we found the maximum α for the extreme values of n (i.e., 0.1 and 0.8); a decrease in α was observed for the intermediate values of n. Figure 8 illustrates examples of good matches between shear-thinning fluid simulations for n = 0.62 and the model based on the fitted shift factor values, given in Tables 8 and 9, for the synthetic silica A1 and Estaillades samples.The results obtained here clearly show that the shift factor increases by increasing the complexity of the porous medium-or, more strictly, the pore space complexity-of the medium.This issue will be discussed in detail later in Section 6.5.
It may be helpful to look at the μ pm -γpm curves for a fixed n vale but varying rock complexity.Figure 9a illustrates simulated porous medium viscosity versus intuitive porous medium shear rate ( γpm = u/√(kϕ)) for different samples but for a fixed value of the flow behavior index (n = 0.62), suggesting that μ pm changes with the pore structure for a given flow or shear rate.This confirms that the shift factor varies with types of porous media (pore-

Table 8
The Estimated Shift Factor (α) From Curve Fitting for Different Values of the Flow Behavior Index, n, for the Power-Law Model Note.The coefficient of determination (R 2 ) for all cases is also reported.space complexity).Figure 9b shows the porous medium viscosity curves for all samples for n = 0.62 when plotted against γpm = αu/√(kϕ), representing that all the viscosity curves of different pore structures coincide with each other.It is evident that the shift factor is a strong function of the pore space complexity and rock type.Our results in Figure 9b are in good agreement with the results published by Tosco et al. (2013) (see Figure 7 of this reference) for the same rheological model (the Cross model with the same parameters) but different porous media.This agreement revealed that the deviations from the decreasing trend in the μ pm -γpm curves observed in this study were due to the onset of the inertial flow effects.

Effect of Flow Behavior Index on the Shift Factor
So far, we have found that the shift factor significantly varies with different pore structures.Figure 11 illustrates the calculated shift factor (α) using the pore-scale flow simulations for different flow behavior indices and various samples with varying pore-space complexity, representing that the larger values of α correspond to extreme Note.The coefficient of determination (R 2 ) for all cases is also reported.

Effects of Rock Pore Complexity on the Shift Factor
Utilizing pore-scale flow simulations provides a unique opportunity to assess the key parameters influencing the complexity of pore space and their impact on the shift factor, which is typically inaccessible through traditional methodologies like laboratory experiments.This section provides the relationship between pore space complexity and the shift factor, utilizing data presented in Figures 12 and 13, which highlight variations for the power-law and the Cross-fluid models.Key observations are as follows: Impact of porosity on shift factor: A clear negative correlation exists between the shift factor and porosity.As the porosity increases, the shift factor significantly reduces.Highly porous rocks exhibit reduced disparity between the bulk and in-situ porous medium viscosities.The underlying reasoning is the probability of a more homogenous viscosity distribution in materials with high porosity.Conversely, rocks with low porosity showcase an elevated shift factor.Such a phenomenon is attributed to the dominance of narrow flow channels and a consequent wide range of shear rates and local viscosities.
Tortuosity's role: A trend emerges when analyzing the shift factor, α, against F.ϕ.The latter serves as a proxy for the tortuosity of pore space.Higher tortuosity corresponds to an increase in the shift factor.This is particularly evident for the Estaillades carbonate sample.Its F.ϕ values entirely surpass those of other samples, as noted in Table 2.In fact, the peak in the shift factor is linked to this sample, underscoring the profound influence of tortuosity.
Cementation exponent: A rising trend in the cementation exponent correlates with an upswing in the shift factor.This hints at the fact that samples exhibiting greater heterogeneity and reduced pore space connectivity show significantly increased shift factors.
Characteristic pore length: Rocks characterized by smaller characteristic pore lengths (analogous to lower permeability) show a larger shift factor.Given this trend, carbonate rocks, often characterized by a substantive microporosity, tend to have shift factors exceeding those observed for sandstones in most cases.
For the power-law fluid flow, the general expressions of the shift factor α based on the capillary bundle model and assuming a constant value for the tortuosity illustrate that α is a function of the fluid rheology (here, the power-law index).However, the literature and also the findings of this work show that α is a function of both the fluid rheology and the pore space complexity.Hence, in line with Equations 12 and 14, it is argued that the coefficient C in the shift factor equation (Equation 15) is not a constant value and is generally a function of the flow behavior index and the pore space complexity.Based on the calculated shift factor using the pore-scale flow simulations, it is possible to estimate the coefficient C. Regarding the Teeuw and Hesselink (T-H) model as the closest model to the simulation data found in Section 4, the calculated C for each sample, from the simulation data, is reported in Table 12.Note that the value of C is equal to √2 ≈ 1.41 from the T-H model for the spherical packs, assuming the capillary bundle model.Balhoff and Thompson (2006) have reported experimental values of C for various non-Newtonian fluids and unconsolidated porous media (packed beds), indicating that it is in the range of 0.8-1.67.It is evident from Table 12 that the values of C for less complex sandstones are close to the results obtained for the unconsolidated porous media, while complex porous media, including carbonate rocks in particular, have significantly higher values of C and hence the shift factor.A close inspection of the effects of fluid rheology and rock types reveals that the value of C increases as the flow behavior index or the complexity of the porous medium increases.Similar results were obtained by Cannella et al. (1988), using an experimental study of Xanthan rheology in a network of capillary tubes with varying tube radii and interconnectivity.In this study, we found that for Sandstone S1 and Estaillades carbonate, the values of C are higher than those predicted from the bundle-ofcapillary-tube models (see Table 1), while synthetic silica A1 shows considerably low values for C and Bentheimer sandstone exhibits an intermediate behavior.

Discussion: Comparison With Literature
Figure 14 shows the pressure gradient against Darcy's velocity, derived from pore-scale flow simulation, and contrasts it with results from the Cross-fluid model, both with and without incorporating the shift factor, for Bentheimer and Estaillades samples.As observed in Figure 14, the Cross-fluid model aligns well with simulated data at lower Darcy velocities, indicating a Newtonian fluid behavior at these velocities.However, as the velocity amplifies, a noticeable disparity emerges between the simulated pressure gradient and the Cross model results that exclude the shift factor.This deviation is particularly pronounced for the Estaillades carbonate sample, Water emphasizing that an accurate shift factor becomes increasingly critical for samples with elevated pore-space complexity.Consequently, the shift factor is indispensable for precisely forecasting the macroscopic pressure differential in non-Newtonian fluid flow within porous rocks.
Table 13 compares the calculated shift factor in this work and those reported in the literature for the Cross-fluid model with exactly the same parameters.Tosco et al. (2013) used 2D SEM images of unpacked sand grains with high porosity and permeability and computed the shift factor.Their results show a relatively similar trend; reduced shift factor for the cases with higher porosity and or permeability, as also evident from the results of Zhang et al. (2019).Nonetheless, we observe some anomalies in data reported by Tosco et al. (2013), which suggests the potential for a non-dimensional analysis scheme with experimental data to develop a relationship that connects the macro and/or microscale properties of the porous media to the shift factor.Furthermore, what is interesting from the data in Table 13 is that the shift factors obtained by Tosco et al. (2013) and Zhang et al. (2019) are in the same range as those calculated in our studies despite wide discrepancies among porous media, representing that the fluid rheology, here the flow behavior index, has a much more significant effect on the shift factor than the porous media properties.

Summary and Conclusions
In this work, we developed a numerical solver based on the OpenFOAM finite volume library to simulate the flow of non-Newtonian fluids in porous media.Our primary objective was to bridge the gap between the behavior of fluids at the macroscopic level and their complex interactions within the complex pore spaces of diverse porous rocks.Our conclusions shed light on various aspects of non-Newtonian fluid flow and its implications for porous  media, contributing to a deeper understanding of this complex phenomenon.While most previous studies have focused on the impact of rheological parameters on flow in porous structures, we have paid particular attention to the influence of pore microgeometry and macroscale rock properties.The main conclusions can be summarized as follows: We found that the shift factor (α), a crucial parameter to accurately predict pressure gradients in non-Newtonian fluid flow in porous media, is significantly influenced by the flow behavior index (n).Interestingly, the maximal shift factor values were observed at extreme n values, while intermediate values of n led to a reduction in α.This dynamic underscores the importance of considering the interplay between flow behavior and pore structure complexity when evaluating the shift factor's impact.At a given porous medium shear rate or flow rate, the porous medium viscosity decreased by increasing the complexity of the pore space.The shift factor decreased by increasing the porosity and the characteristic pore size.However, it increased by increasing pore tortuosity and heterogeneity and decreasing pore connectivity.Despite the simplicity of the capillary bundle model, this approximation can correlate the viscosity of non-Newtonian fluids in porous materials, provided that proper coefficient values are used, for example, from pore-scale flow simulation.
Finally, we remark on the applications and limitations of the methods and results achieved by this work.The method adopted in this study and, in general, the pore-scale simulation approach is a strong, versatile alternative to laboratory tests and analytical methods.This technique allows exploring various fluid models, a wide range of flow rates and pressures, and different 3D microstructures.This study showed that the pore-scale simulation approach is useful in accurately predicting pressure drop in the flow of non-Newtonian fluids through porous materials by calculating an accurate shift factor for a given fluid-porous medium system.However, the effectiveness of this method depends on knowing detailed information about the pore space of the medium, which can be obtained through high-resolution 3D microtomography imaging.This limitation restricts the applicability of the pore-scale numerical simulation to μ-CT images of complex tight porous media.The abundance of micro-and nanopores combined with pore structure complexity is a critical issue that needs further investigation.Another issue of worthy further investigation is evaluating the computational performance of the proposed solver for simulating non-Newtonian flow in large data sets where massive parallel computing is required.

Table 13
A Comparison Between the Calculated Shift Factor in This Work and Those Obtained by Tosco et al. (2013) and Zhang et al. (2019)   Water Resources Research 10.1029/2023WR036125 Newtonian and non-Newtonian flow simulations were applied on 3D voxelized pore space images of four samples obtained from the open-access μ-CT image library available at the Imperial College Water Resources Research 10.1029/2023WR036125 London website (https://www.imperial.ac.uk/earth-science/research/research-groups/pore-scale-modelling/micro-ct-images-and-networks/).The images used in this study are (a) a synthetic silica sample named A1, (b) the Bentheimer sandstone, (c) a sandstone named S1, and (d) the Estaillades carbonate (Figure 1).Table
in terms of μ pm against the intuitive porous medium shear rate, defined by γpm = u/√(kϕ) for different flow behavior indices (n = 0.2, 0.4, 0.6, and 0.8).The calculated shift factor for each n is also shown in Figure2.The mean shift factor for the power-law flow with varying flow behavior index through the beadpack sample is 1.57 ± 0.03.It is worth mentioning that the power-law equation represents a linear relation in the log-log plot of the porous medium viscosity versus shear rate.The deviation from the straight line at high flow rates (corresponding to high shear rates) observed in the simulation data of Figure2marks the start of the inertial flow effects.For low values of n (i.e., for the cases of n = 0.2 and 0.4), numerical results show a deviation from the linear relationship (creeping flow region) at extremely low flow rates, probably due to the importance of the interaction between the fluid and the pore walls, which violates the assumptions of the linear Darcy's law model reveals that the numerically simulated μ pm varies quite a bit and is in between the results of the C-M (Christopher and Middleman) and T-H (Teeuw and Hesselink) models, except for n = 0.2 where the simulation data is between C-M and Savins models.It can be argued that as the power-law index increases, the numerical results are closer and match the T-H model.It is interesting to check the numerical solver

Figure 2 .
Figure 2. Comparison of semi-analytical models (Equation14) of the power-law fluid flow in a beadpack with our pore-scale numerical simulation results in terms of porous medium viscosity, μ pm , as a function of the intuitive porous medium shear rate, defined by u/√(kϕ).Note that the value of K (the flow consistency index) was taken to be 0.3 Pa.s n .The values of ϕ and k are reported in Table2.

Figure 3 .
Figure 3.Comparison of the calculated α values from pore-scale flow simulations for different values of n with those obtained for semi-analytical models (Equation14and Table1).

Figure 4 .
Figure 4.A comparison between the experimental and simulated values of α/√(kϕ); The experimental data of (a) and (b) were obtained from Chauveteau (1982) and Rodríguez de Castro and Radilla (2017b), respectively.The solid gray line indicates the equidistant line from both the x-and y-axes.The dashed lines illustrate the prediction range at the confidence level of 95%.

Figure 5
Figure 5 illustrates the dimensionless pressure drop ratio (PDR) for Bentheimer and Estaillades as functions of Reynolds numbers (Equations 18 and 23) for the flow of Newtonian fluid and non-Newtonian fluids (n = 0.62 and 0.8) depicting the transition from creeping to inertial flow regimes.Muljadi et al. (2016) adopted the transition point for the creeping flow regime as the point at which the pressure drop resulting from the linear Darcy's law fell

Figure 5 .
Figure 5.The ratio of pressure drop at creeping flow to the apparent pressure drop as a function of Reynolds number depicting the transition from the creeping to inertial flow regimes for Bentheimer and Estaillades samples for flow of Newtonian, power-law, and Cross fluids.

Figure 6 .
Figure 6.The results of pore-scale flow simulations in terms of porous medium viscosity (μ pm = k/u • Δp/L) versus the intuitive porous medium shear rate ( γpm = u/ √(kϕ)), for the power-law rheological model for different values of the flow behavior index.

Figure 7 .
Figure 7.The results of pore-scale flow simulations in terms of porous medium viscosity (μ pm = k/u • Δp/L) versus the intuitive porous medium shear rate ( γpm = u/ √(kϕ)), for the Cross fluid rheological model for different values of the flow behavior.

Figure 10
Figure 10 illustrates the effects of the time constant of the Cross-fluid model (m) on the porous medium viscosity versus the intuitive porous medium shear rate.Similar behavior was found for all samples where μ pm decreased with increasing the time constant at all shear rates.Table 10 summarizes the calculated values of the shift factor for different values of the time constant, m, for the Cross-fluid model.Table 11 summarizes major statistics (the mean and coefficient of variation) for the simulated shift factor values for different values of flow behavior index, n, and the time constant of the Cross model, m.The coefficient of variations data indicates that the variability of the shift factor is higher by changing the flow behavior index compared to the variability due to the time constant.Hence, it can be argued that the flow behavior index is a more important rheological parameter affecting the shift factor.

Figure 8 .
Figure 8.The simulation data of μ pm versus γpm and the model based on the fitted shift factor, given in Tables8 and 9, for the synthetic silica A1 and Estaillades for the power-law and Cross-fluid models with n = 0.62.

Figure 9 .
Figure 9. Variations of porous medium viscosity (μ pm ) with porous medium shear rate ( γpm ) (a) without and (b) with shift factor for samples with various complexities.The pore-scale numerical simulation data are for the Cross-fluid model with n = 0.62.

Figure 10 .
Figure 10.The results of pore-scale flow simulations in terms of porous medium viscosity (μ pm = k/u • Δp/L) versus the intuitive porous medium shear rate ( γpm = u/ √(kϕ)), for the Cross fluid rheological model for different values of the time constant.

Figure 12 .
Figure 12.Variations in average shift factor as functions of porosity, the product of formation factor and porosity (as an indication of pore tortuosity), the cementation exponent, and the characteristic pore length for the power-law fluid model.

Figure 13 .
Figure 13.Variations in average shift factor as functions of porosity, the product of formation factor and porosity (as an indication of pore tortuosity), the cementation exponent, and the characteristic pore length for the Cross-fluid model.

Figure 14 .
Figure 14.Plots of pressure gradient against Darcy's velocity for Bentheimer sandstone and Estaillades carbonate samples.The pore-scale numerical simulation data are for the Cross-fluid model with n = 0.62.
for the Cross-Fluid Model With the Same Parameters ( area (based on the volume of the porous medium).

Table 1
The Values of Constants C and D in Equation 14 and the Corresponding Formulas of Calculating Shift Factor α for the Power-Law Flow Through Porous Media Based on the Capillary Bundle Model

Table 2
The Image Properties and the Calculated Parameters From Single-Phase Newtonian Flow Simulations a computed by √k.b m c = log (F)/log (ϕ).

Table 2 .
It is interesting to mention that from Table 2, a perfect agreement is observed

Table 3
Computational Resources Used in This Study and Simulation Times AMIRI ET AL.

Table 4
The Values of Parameters in the Carreau Model Used in the Simulations for Comparison With the Experimental Data

Table 6
Values of Parameters in the Power-Law and Cross-Fluid Models of Shear-Thinning Fluids Used in Our Simulations AMIRI ET AL.

Table 7
The Estimated Critical Reynolds Number (Onsets of the Inertial Flow Regime) for the Samples Used in This Study for the Cases Where n = 0.62

Table 9
The Estimated Shift Factor (α) From Curve Fitting for Different Values of the Flow Behavior Index, n, for the Cross-Fluid Model

Table 10
The Estimated Shift Factor (α) From Curve Fitting for Different Values of the Time Constant, m, for the Cross-Fluid Model Note.The coefficient of determination (R 2 ) for all cases is also reported.

Table 11
The Arithmetic Means and Coefficients of Variation of the Calculated Shift Factor for Different Power-Law Indices (for Both the Power-Law and Cross-Fluid Models) and for Different Time Constants of the Cross-Fluid Model Figure 11.Values of shift factor for different samples and flow behavior indices.

Table 12
The Calculated Coefficient C in the Shift Factor Equation (Equation15) Based on the Results Obtained by the Pore-Scale Flow Simulations for the Power-Law Fluid Flow and theTeeuw and Hesselink Model