Upscaling Hydrodynamic Dispersion in Non‐Newtonian Fluid Flow Through Porous Media

Hydrodynamic dispersion in flow through porous media is an essential phenomenon in many geosystems. While dispersion in flow of Newtonian fluids is relatively well‐understood, many subsurface applications, such as groundwater remediation, use the flooding system with non‐Newtonian fluids, such as polymer oxidants. Despite its significance, however, very limited studies have been carried out focusing on hydrodynamic dispersion in flow of non‐Newtonian fluids. The present study addresses a fundamental question regarding how solute transport in flow of non‐Newtonian fluids in porous media differs from the Newtonian limit. We report on the results of an extensive pore‐scale study of the upscaling of advection‐diffusion in flow of a non‐Newtonian, shear thinning fluid through disordered and spatially‐correlated porous media, using an advanced GPU‐based pore‐scale simulator in order to delineate the effects of the fluid's rheology, dynamics of fluid flow, and pore‐scale spatial correlations on hydrodynamic dispersion. The simulations indicate a surprising non‐monotonic relation between the injection rate (or injection velocity) and the dispersivity. While dispersivity in Newtonian fluid flow in porous media is constant in the flow regime that we study, we find, however, that the dispersivity in flow of non‐Newtonian fluids is shear‐dependent. This highlights the gap in the existing theories of transport through porous materials for non‐Newtonian fluids, which can lead to erroneous estimates of dispersion coefficients in porous materials when the fluid's rheology cannot be represented by Newtonian mechanics.

Non-Newtonian fluids exhibit complex rheological properties (Bird et al., 1987), such as shear-dependent viscosity, viscoelastic behavior, and time-dependent rheological characteristics, usually referred to as thixotropy and antithixotropy. A non-Newtonian fluid is shear thinning (ST) or pseudoplastic, if its apparent viscosity μ decreases with increasing shear stress τ, and it is shear thickening or dilatant when μ increases with increasing τ (Bird et al., 1987). The former type is quite common and includes blood, ketchup, syrup, and crude oil (Sochi, 2007(Sochi, , 2010. In fact, it was recently shown that even water in small tubes or pores could exhibit ST behavior (Cobeña-Reyes & Sahimi, 2020). To describe the rheology of non-Newtonian fluids, several theoretical, empirical, and semiempirical models have been proposed (Bird et al., 1987;Shende et al., 2020).
The interplay between the complex rheology of non-Newtonian fluids and the disordered morphology of porous media implies that the upscaling of flow and transport of non-Newtonian fluids in such media is much more complicated than that for Newtonian fluids (Sahimi, 1993;Zami-Pierre et al., 2017). It is, for example, clear that, due to the significant spatial variations in the pore sizes and pore-scale flow velocities , the pore-scale shear stresses and, hence, the pore-scale viscosity of a non-Newtonian fluid also varies spatially. Thus, the rheological characteristics of bulk non-Newtonian fluids are not directly applicable to the flow of the same fluids in porous media. Chiapponi et al. (2020) experimentally and numerically studied the flow and transport in a Hele-Shaw cell and concluded that the fluid dispersion is very sensitive to the rheology of non-Newtonian fluids even in a layered fracture. Therefore, the critical questions of how the properties of non-Newtonian fluids in porous media should be upscaled, and how the representative elementary volume (REV) can be identified, must be addressed. An effective viscosity is usually introduced in order to account for the nonlinear effect of the fluid's rheological characteristics and the flow's heterogeneity (Eberhard et al., 2019(Eberhard et al., , 2020Seybold et al., 2021). Several analytical, empirical, and semiempirical correlations have been developed for idealized models of porous media in order to determine a correction factor for estimating the effective viscosity in porous media (Delshad et al., 2008) from data for bulk viscosity. Depending on the fluid and the morphology of the pore space, however, such correction factors exhibit up to three orders of magnitude variations, which is not supported by experimental data (Berg & van Wunnik, 2017).
Study of the ST fluids in porous media entails, as the first step, understanding their flow in straight capillaries (Metzner, 1957). Using a generalized Darcy's law, Sadowski and Bird (1965) determined an apparent viscosity for shear-sensitive fluids in tubes. Dependence of the pressure drop on the flow rates in inertial ST flow in a packing of beads was investigated experimentally and by numerical simulation (de Castro & Radilla, 2017;Sabiri & Comiti, 1995). Sabiri and Comiti (1995) used computational simulation to study the effect of tortuosity and pore sizes on the pressure drop in the flow of ST fluids through a packed bed. Morais et al. (2009) simulated the flow of shear thinning fluids in 3D disordered porous media with large porosity, and proposed a universal curve for describing the relation between the effective permeability and a suitably-defined Reynolds number. Flow of a ST fluid through an array of cylinders was studied by Zami-Pierre et al. (2018), who reported that disorder reduces the effect of the nonlinearity on the average flow. Sahimi (1993) used the critical path analysis to derive an expression for the effective permeability of power law fluids, including ST ones, in a porous medium.
Pore-network (PN) models have been used widely to study the flow of non-Newtonian fluids, particularly ST ones (Lopez et al., 2003;Sorbie et al., 1989), in porous media. Balhoff and Thompson (2004) studied the flow of guar gum solution in a PN model of packed beds, while Sochi (2007) simulated single-phase flow of three types of non-Newtonian fluids and stated that the ST accentuates the effect of the heterogeneity. Using a 3D PN, Xiong et al. (2020) addressed the impact of fluid rheology and pore connectivity on the permeability, while Didari et al. (2020) studied the flow of a more complex non-Newtonian fluid, the Bingham model, in a 3D PN model. Pearson and Tardy (2002) reviewed models of single and multiphase transport in non-Newtonian fluid through porous media and concluded that PN models provide the most reliable estimates of various properties. Tracer transport in Newtonian fluids was also studied using PN models (Mehmani et al., 2015;Mehmani & Tchelepi, 2017;Yang et al., 2016), with some considering the stream splitting in order to account for pore-scale mixing in pore throats (Mehmani et al., 2014). In addition, recently, An et al. (2020) studied solute transport in unsaturated porous media by coupling PN models with a GPU-accelerated numerical algorithm to enlarge the simulation domain and higher computational speed, which made it possible to study the effect of heterogeneity at the REV scale.

Validation: Senyou An, Vahid Niasar
Writing -original draft: Senyou An Writing -review & editing: Senyou An, Muhammad Sahimi, Vahid Niasar But, despite its significance, solute transport in flow of non-Newtonian fluids in porous media has not been addressed adequately. The present study aims to highlight the key difference between the hydrodynamic dispersion in flow of non-Newtonian fluids in porous media and their Newtonian counterparts, and the role of spatial heterogeneity of porous media in the phenomenon. In this paper, we use the aforementioned GPU-based PN model to simulate the flow and transport of a ST fluid in a porous medium. Our primary goal is not the simulation of the phenomenon in a specific porous medium, for example, the Berea sandstone, since we are not aware of any experimental data to compare with the predictions of our model, rather our goal is to lay the foundation for a rigorous study of the physics of the phenomenon. As a result, the methodology that we describe is completely general and applicable to the dispersion in flow of many other types of non-Newtonian fluids in porous media. We use a modified Meter model (Meter & Bird, 1964) (MMM) to describe the fluid's rheology and simulate spatially-correlated or uncorrelated PNs over length scales much larger than the REV, on the order of several centimeters, in order to study the upscaled viscosity and the dispersion coefficient for a ST fluid.
We point out that although simulation of fluid flow and transport in 2D or 3D images of porous media, instead of PN models, is becoming quite common, such studies (e.g., lattice-Boltzmann (LB) based approaches (Boek et al., 2003;Sullivan et al., 2006)) are developed for flow simulation but not yet well-developed for simulation of solute/heat transport in non-Newtonian fluids, usually with expensive computational cost. Naturally, if experimental data for a specific porous medium, plus its image were available, and if the LB method could be developed for a class of non-Newtonian fluid, we could carry out the same type of study with the LB, or extract the PN model of the pore space from its image. But, in the absence of such data and the associated image of the pore space, the methodology that we develop in this paper offers a realistic approach to study the problem of interest.
The rest of this paper is organized as follows. In the next section, the methodology and governing equations for the PN model and upscaling the flow and transport are explained. We then present the results for various injection rates and spatial correlation lengths for a given shear thinning fluid. Finally, we discuss the implications of upscaling flow and transport of non-Newtonian fluids in porous media.

Network Generation
Generation of the spatially-correlated PNs is described in Supporting Information S1. Briefly, an identical network topology was also developed by randomly generating the locations of the pore bodies, which were connected by the Delaunay triangulation method. Then, the sizes of pore bodies and pore throats were selected from distinct correlated fields. The correlated fields were generated using the fast Fourier transform method in order to obtain the spectral density functions, given an autocorrelated field with anisotropic exponential variograms. The pore throats were assumed to be cylindrical, whereas the pore bodies were spherical.
Three-dimensional unstructured PNs in three different shapes were generated for different objectives (see below). A cylindrical PN was generated to validate the model against bead-packing experiments in Section 3.1. The majority of the results are obtained for a radial flow system in a thin disk PN to visualize the front movement better. The PN contained about one million pore bodies and four million pore throats, with an average pore radius of 2.57 × 10 −5 m and a variance of 2.12 × 10 −10 m 2 . The radius R of the PNs is 50 mm and its thickness is 0.2 mm. Based on the average size of the pore body and the size of geometry, the thickness and radius of the sample are roughly equal to 8 and 2000 times of pore radius, respectively. We generated both uncorrelated and correlated PNs with correlation lengths λ/R = 0, 1/40, and 1/20; see Figures 1a-1g. Moreover, to ensure that the results do not depend on the flow pattern, a linear flow system was simulated in a cuboid PN, as shown in the Supporting Information S1.

Pore-Network Model of Flow and Transport of a Non-Newtonian Fluid
We computed the pressure field for the flow of a ST fluid in the PNs. Mass balance at each pore unit-a pore body and half of the connected pore throats-was satisfied by solving the continuity equation, Equation 1a below. We assumed the Hagen-Poiseuille flow, Equation 1b, to be valid, except that the effective viscosity varies from pore to pore, which has been shown to mimic the well slow flow of non-Newtonian fluids in tubes. Note that the throat conductance k ij depends on the effective viscosity, which we refer to the viscosity averaged over a single-pore throat. The apparent viscosity is used to name the upscaled average viscosity in a porous medium. Equation 1a, when written for all pore bodies, must be solved iteratively to reach convergence.
where r ij and l ij are the radius and length of throat ij, and ρ ij are the effective viscosity and density of the fluid in the throat, ΔP ij is the pressure difference between the two ends of the throat, g is gravitational acceleration, z i is the absolute height of pore body i, and V i is its volume. The effective viscosity of the ST fluid depends on the shear stress. Experimental data (the bulk rheology data used in this work were adopted from Xue and Sethi (2012), which had been obtained using the Couette-flow geometry) for the entire range of the viscosity variations were fitted to the modified Meter model (MMM), Equation 2, by minimizing the relative error.
where μ(τ) is the shear stress-dependent viscosity of the ST fluid at shear stress τ, μ 0 , and μ ∞ are the viscosities at zero and infinite shear stress, τ m is shear stress of the ST fluid at which μ drops to (μ 0 + μ ∞ )/2, and S τ is a shear stress-dependent exponent.
The constitutive relation for many non-Newtonian fluids can be accurately modeled with the Meter model, as described by Equation 2. In cylindrical throats of the PN, the shear stress, τ = βr ij ΔP ij /2l ij , is incorporated into Equation 2 in order to calculate the effective viscosity of a ST fluid.
(3) Here, β is a geometrical correction factor to relate the velocity profile in a non-Newtonian fluid using the MMM parameter. Given that for a non-Newtonian fluid, viscosity changes across a pore, β is a geometrical factor for upscaling the viscosity for a given pore. Following Shende et al. (2021) β = 0.8 was used for circular capillary tubes. Solving Equations 1b and 3 iteratively, we obtain the steady state velocity field of the ST fluid in the PN, which is then used to describe the transport of a neutral solute in the PN using the advection-diffusion equation, Equation 4a: where c i is the solute concentration in pore unit i, q ij is the volumetric flow rate in a throat ij, N i is the coordination number of pore body i, eff is the effective dispersion coefficient in a capillary tube for a power law fluid (Dejam, 2018;Vartuli et al., 1995), and D m is the molecular diffusion coefficient. Here, Pe ij = 2v ij r ij /D m is the pore-scale Péclet number, with v ij being the average flow velocity in the pore throat ij, and n is the exponent of the power law, with 0 < n < 1 for a ST fluid, and n = 1 for Newtonian fluids. The data for the linear part of the viscosity-shear rate relation in logarithmic form was fitted using the power law relation, μ(γ) = aγ n−1 , with γ being the shear rate, and a a fitted constant. The modified Meter model that we used is a combination of a power law model with lower and upper cutoffs. During the simulation, if the viscosity of throat belonged to the power law part, the fitted n 0 was used to modify the dispersion coefficient. Whereas in the throats, in which viscosity was in the transition part from power law to plateau section, n was interpolated between n 0 and one based on the slope, as shown in Figure 2b.
Equations 1a and 1b, when written for all the pore bodies, give rise to a set of equations that we solve by a linear matrix solver using an in-house Jacobi preconditioned biconjugate gradient solver, which is fully parallelized by CUDA. Simultaneously, Equation 3 is explicitly updated during the flow simulation. Equations 4a and 4b are explicitly solved for all the pore bodies at every time step based on the velocity field and concentration distribution computed in the last time step. Time step for each pore i is determined based on the residence time, the time required to reach a concentration of 1 for loading or 0 for unloading pores. The global time step is the smallest local time step, considering both loading (increasing concentration) and unloading (decreasing concentration) in pore units.

Boundary Conditions
The simulations in all the geometries were carried out under a constant flow rate at the inlet by adding one more equation to the linear system of Equation 1b, given by where Q is the injection flow rate, m is the number of throats contacted to the inlet reservoir, P 0 is the transient reservoir pressure, P i is the pressure of the inlet pore body i, and k i0 is the flow conductance of the inlet throat. At the outlet, the pressure was set to zero. We carried out the simulations with Q = 10, 1.0, 0.1, and 0.001 mm 3 /s.
For the simulation of solute transport, the concentration of the inlet reservoir was assumed to be one. For the outlet boundary pores, the fully-developed boundary condition was utilized, implying that there was no diffusion contribution to Equation 4a between outlet pores and outlet reservoir (∂C/∂x = 0). The molecular diffusion coefficient D m was 1.0 × 10 −9 m 2 /s.

Computational Setup and GPU-Based Acceleration
The thin disk networks with specific correlation lengths, as well as the given pore size distribution, were generated based on the method introduced in Section 2.1. As aforementioned, the properties of the non-Newtonian fluid and the governing equations for flow and solute transport were described. To validate the proposed algorithm, a cylindrical network was extracted to represent a randomly packed bed of monodisperse spherical beads, and the boundary conditions imposed on the PN were kept consistent with the experimental setup.
To accelerate the computational speed and enlarge the simulation domain, flow and solute transport simulations were fully parallelized using the proposed GPU-CUDA algorithm (An, Yu, Wang, et al., 2017;An et al., 2022). A GPU-based linear solver based on Jacobi's preconditioned conjugation gradient method was adopted to solve the mass balance equations . The governing equations for solute transport through the porous medium were explicitly solved by allocating pore bodies to various kernels to expand the loop calculation. Both the flow and transport simulations were carried out on the NVIDIA Tesla v100 GPU card, which has 5120 CUDA cores with 1,380 MHz clock frequency, 900 GB/sec memory bandwidth, and 32 GB of global memory.

Modeling Parameters
We fitted the experimental data for Xanthan gum with a concentration of 6.0 g/L in water at 25°C (Xue & Sethi, 2012) to the modified Meter model to obtain a relationship between the viscosity and shear stress; we obtained, μ 0 = 100.6 Pa⋅s, μ ∞ = 6.2 × 10 −3 Pa⋅s, τ m = 0.75 Pa, and S τ = 3.4, shown in Figure 2a. As mentioned earlier, in Equation 4b n is the exponent of the power law model, which for our data is n 0 ≃ 0.232.

Estimation of Upscaled Parameters
Given the radial geometry of the simulated system, the absolute permeability of the PN was calculated as follows: where μ e is the effective viscosity of the fluid in the pore space, k is the absolute permeability, h is the thickness of the PN, ΔP is the pressure difference between the inlet and outlet, and R i and R o are, respectively, the radii of the inlet and outlet. By comparing the simulation results for the non-Newtonian fluid with those for the Newtonian one, the effective viscosity of the former fluid in the PN is estimated.
To quantitatively analyze the concentration distribution, the effective dispersion coefficient is fitted using an analytical equation for the solution of concentration distribution as a function of time in a radial flow, which was derived by Lau (1959); see also Sauty (1980). If the solute is injected continuously into the PN at its center, and a fully-developed condition is assumed at the outlet, the approximate solution for the dimensionless concentration C R is given by where u(R) is the average pore velocity at a distance R from the center, which is the radius of the disk, D is the effective dispersion coefficient, and τ R is a dimensionless time, τ R = Qt/ϕV, in which V is the volume of the PN, and ϕ is the porosity.
To estimate the dispersivity α, we used D = αu β . We assumed that β = 1, and estimated α for various limits, using the average pore velocity at the outlet of the domain u . In the present work, the PN results were fitted to Equation 7 in order to calculate the effective dispersion coefficient and dispersivity.
To quantitatively describe the concentration patterns in the pore space, we define a shape factor S, S = A/P 2 , where A is the area and P is the perimeter of the patterns. Similar analyses were performed for linear flooding, as presented in the Supporting Information S1.

Results and Discussion
As mentioned earlier, we carried out PN simulations in both correlated and uncorrelated PNs for a Newtonian (water) and a ST fluid (Xanthan gum). The absolute permeabilities of the uncorrelated, moderately-correlated, and highly-correlated PNs were 0.086, 0.117, and 0.125 Darcy, respectively. In what follows, we first present the results for validating the results for the ST fluid flow, and then analyze the spatial and temporal distribution of various quantities in order to understand the effect of flow dynamics and pore-size correlation length on solute transport in flow of both Newtonian and non-Newtonian fluids through porous media.

Validation by Comparing With Experimental Data
To validate whether the proposed model can accurately simulate non-Newtonian fluid flow in porous media, we compared the simulation results of our model with experimental data reported in Eberhard et al. (2019). The porous medium was made of packed, monodisperse spherical beads.
The network was extracted from the generated porous packed bed, shown in Figure 3a. To accurately fit the absolute permeability of the bead packing using the pore-network model, the pore sizes were increased by a factor of 1.07. The measured bulk rheology was fitted using the MMM as shown in Figure 3b. With the bulk rheology data and pore-network morphology, the flow of the non-Newtonian fluid in the PN was simulated under isothermal conditions and with various injection rates. As a result, we established the relation between the Darcy velocity and the apparent viscosity of the non-Newtonian fluid in the porous medium and compared the results with the experimental data as shown in Figure 3c. The results clearly demonstrate that the model is capable of reproducing experimental data.

Uncorrelated Network
Flow and solute transport were simulated in an uncorrelated PN. The results for average concentrations C R of 0.4 and 0.8 and an injecting flow rate of 1.0 mm 3 /s are shown in Figure 4. The concentration profiles in both Newtonian and ST fluids have smooth fronts, as the medium is uncorrelated. Furthermore, for each specific flow rate, the solute transport in Xanthan gum requires a longer time than in water to reach a given average value after the breakthrough point, since the flow of the ST fluid generates slight fingering in the concentration pattern.
We fitted the concentration profiles to Equation 7 in order to estimate the dispersivity, as shown in Figure 4d. The fitting accuracy is shown in the Supporting Information S1. Due to the fitting error, the fitted dispersivity length for Newtonian (water) flow has a slight fluctuation. The dispersivity values were averaged to reflect the fact that the dispersivity is a property of the porous medium when saturated by Newtonian fluid in the linear flow region and does not depend on the flow rate. The dispersivity varies only weakly with increased flow rates, due to the uncorrelated structure of the PN and similar velocity distribution. Moreover, the distribution of the pore shear stresses, τ ij = βr ij ΔP ij /2l ij , was utilized for analyzing solute transport. For a flow rate of 10.0 mm 3 /s, the viscosity in most of the pore throats approached μ ∞ , with the apparent viscosity being 6.7 mPa⋅s (see Table 1), and the corresponding average shear stress being 43.0 Pa, close to the average value shown in Figure 4c. Decreasing the injection flow rate led to a smaller shear stress, whereas the apparent viscosity was much larger than its average value. This is due to the fact that the fluid in pore throats that experience high shear stress would have a smaller viscosity and dominate the flow pathways in the case of the ST fluid.

Correlated Networks
Flow and solute transport were simulated in spatially-correlated PNs with correlation lengths λ = 1.25 and 2.5 mm and a network radius of 50 mm. In the moderately-correlated PN with λ = 1.25 mm, fluid flow was first simulated to obtain the velocity and shear stress distributions, as shown in Figure 5. The distributions of local shear stress in pore throats are shown in Figure 6b. Similar to flow in the uncorrelated PN, the apparent viscosity is close to its average with a flow rate of 10.0 mm 3 /s, as the viscosity in most pore throats is close to μ ∞ . For other flow rates, the shear stress at the apparent viscosity is much larger than the average shear stress in the entire domain, given that not all the pore space contributes equally to the total flow. The apparent viscosity was 6.9, 11.7, 41.3, and 1377.1 mPa⋅s for the flow rate of 10, 1.0, 0.1, and 0.001 mm 3 /s, respectively, with the corresponding bulk shear stress being 31.2, 13.7, 7.8, and 2.6 Pa. Compared to the average shear stress in Table 1, the bulk shear stress has a larger difference from its average value when the flow rate decreases.
Having computed the steady-state velocity field, solute transport was simulated in the PN in order to determine how the concentration distribution evolves with time; see Figure 6a. When the Newtonian fluid is injected into the pore space, the concentration distributions exhibited slight differences for the four flow rates, as seen in Figure 6a. The concentration patterns differ significantly from those in the Newtonian fluid for a PN saturated by a ST fluid. For a given concentration, the required time to reach that concentration exhibits a non-monotonic   Note. Q is the flow rate, v is the average pore velocity, is the average shear stress, μ B is the bulk viscosity at , and represents the average viscosity. To quantitatively demonstrate the non-monotonic relation, we calculated the dispersivity for various flow rates. As shown in Figure 6d, the dispersivity reaches a maximum value for Q = 0.1 mm 3 /s. To further analyze the  concentration patterns, the shape factor of the concentration fronts at an average concentration of 0.4 was calculated; see Figures 5 and 6c. They indicate that the concentration pattern manifests the most severe fingering, Q = 0.1 mm 3 /s, since diffusion, a stochastic process, is dominant.
Flow and solute transport were also simulated in the highly-correlated PN with a correlation length, λ = 2.5 mm; see Figure 7. In the case of the Newtonian fluid, the increase in the injection flow rate resulted in only a slightly sharper front, whereas for the ST fluid, the non-monotonic behavior was further enhanced with the λ. The disper sivity grows from 25.8 to 37.6 mm when the flow rate Q increases from 0.01 to 0.1 mm 3 /s and decreases to 18.9 mm when Q = 10.0 mm 3 /s. The shape factors of the concentration fronts at the average concentration of 0.4 also exhibit a non-monotonic behavior. The simulations indicate that solute transport in flow of a ST fluid through the PNs produces non-monotonic dependence of the dispersivity on the flow rate, and more complex fingering patterns, with the latter quantified with a flow-rate dependent shape factor. The non-monotonic behavior is due to the fact that at very small (high) injection rates, the pore-scale shear stress and shear rate are mostly associated with the low (high)-range of the plateau regions of the bulk rheology, where viscosity is not shear-dependent anymore. Thus, at low or high injection rates, there is practically no spatial variation of the viscosity, and the ST fluid behaves similarly to a Newtonian fluid. Thus, dispersivity should follow almost what is obtained for a Newtonian fluid. However, for intermediate injection rates, pore-scale shear-stress/shear-rate encompass a broad range of viscosity variations, from μ min to μ max . As a result, pore-scale viscosity varies spatially across the PN. For a ST fluid, the viscosity decreases, hence the pore-scale velocity and the dispersion coefficient both increase. If we define the dispersivity α with D = αv β , it increases with increasing v if β = 1. However, assuming α is a length scale and a network property in the studied regime, we conclude that β should change for ST fluids when v varies.
The interplay between the bulk rheology, the dynamics, the injection rate, and the porous medium heterogeneity would determine the spatial pattern of solute transport and fingering, as well as the apparent viscosity. If the porous medium is anisotropic, the anisotropy will also play a role.

Conclusions
This paper addressed how solute transport in flow of ST fluids through porous media differs from that with Newtonian fluids. Three-dimensional unstructured PNs with three different shapes were generated based on specific statistical information for the pore sizes and the correlation lengths. A cylindrical PN was used to validate the model against bead packing experiments. The majority of the results were obtained for a radial system in a disk-shaped PN to better visualize the front movement. Similar results were obtained for a linear flow system in a cuboid PN, as shown in the Supporting Information S1.
Single-phase flow of non-Newtonian fluid and solute transport in the flow field were then simulated. To include both the power law and plateau regions in the viscosity diagram, the experimental data for the dependence of the viscosity on the shear stress were fitted to the modified Meter model. With the GPU-accelerated algorithm, we simulated solute transport in large PNs and included spatial heterogeneity in the PNs that had a radius of 0.5 cm with one million pore bodies.
Solute transport in flow of a ST fluid through the PNs produces non-monotonic behavior for dispersivity and fingering patterns, with the latter quantified with a shape factor as a function of the flow rate. The dispersivity length for Newtonian (water in this case) flow is a property of the porous medium and does not depend on the flow rate in the specific regime. For the ST fluid, the non-monotonous dependence of the dispersivity length on the imposed flow rate is due to the fact that at very small and very large flow rates the flow behaves similar to Newtonian fluid. Thus, dispersivity should almost follow what is obtained for a Newtonian fluid. For intermediate injection rates, however, pore-scale shear-stress/shear-rate encompass a broad range of viscosity variations, from μ min to μ max . As a result, pore-scale viscosity varies spatially across the PN. For a ST fluid, the viscosity decreases, hence the pore-scale velocity and the dispersion coefficient both increase.
Our results demonstrate that the spatial correlations enhance fingering in ST fluids. Clustering of the high-and low-permeability regions allows the development of preferential flow paths, such that most of a ST fluid is transported through the high-permeability regions. As a result, the bulk rheology of the fluid cannot be directly utilized for estimating the upscaled porous media viscosity. It is essential to capture the shear distribution in the preferential flow paths in order to be able to estimate the upscaled viscosity with reasonable accuracy.

Data Availability Statement
The readers can find the input data and raw data of figures presented in the manuscript in the Mendeley repository with the An (2022, doi: 10.17632/tczrhyrsj6.1).