Polymer extension at stagnation points governs flow thickening of polymer solutions in ordered porous media

The flow of viscoelastic polymer solutions through porous media at low Reynolds number ($\operatorname{Re} \ll1$) underlies diverse chemical production [1]–[4], additive manufacturing [5, 6], separation [7]–[10], energy production [11]–[18], and environmental remediation [19]–[22] processes. And yet, predicting how pore-scale dynamics set the macroscopic flow resistance, quantified by the “apparent viscosity” $\eta_{\mathrm{app}} \equiv\frac{\Delta P}{L}\frac{k}{Q/A}$, remains an open problem; here, Darcy’s law relates the total pressure drop $\Delta P$ across a length $L$ of the medium, the volumetric flow rate $Q$ through the cross section $A$, and the permeability $k$ [23]. Indeed, despite the typical shear-thinning rheology of these fluids, they anomalously ‘flow thicken’ in a porous medium: $\eta_{\mathrm{app}}$ increases, not decreases, above a threshold flow rate [24]–[32]. The underlying reason has been a mystery for over half a century. Various mechanisms have been proposed [2, 11, 13, 29]–[45], [47]–[?], [?]–[?], [?]–[?], [?]–[?], [?]–[?], but quantifying their relative contributions to flow thickening has been challenging due to the opacity of typical media. Instead, prior work has relied on numerical simulations to characterize pore-scale flow dynamics [51], [52], [?], [?], [?]–[?] or on semi-empirical models with adjustable fitting parameters [?]–[?].

Recent work has identified the central roles of two key mechanisms. Direct flow visualization within disordered 3D sphere packings [2, 48] demonstrated that the onset of an elastic flow instability—which generates unsteady pore-scale flow fluctuations [6], [11], [43], [?], [?]–[?], [?]–[?], [?]–[?], [?, ?]—coincides with the onset of flow thickening, and that the resulting enhancement of viscous dissipation accounts for much of the increase in $\eta_{\mathrm{app}}$ near this onset. However, the elastic instability contribution alone progressively underpredicts flow thickening at higher flow rates—an effect that is particularly pronounced in ordered 3D sphere packings, where polymeric stresses generated by extensional flow have been suggested to govern flow thickening instead [40], as corroborated by independent 2D experiments [?]. Yet a quantitative way to link extension at the pore scale to macroscopic flow thickening in porous media has remained lacking. Here, we develop this link from a mechanical energy balance and validate it experimentally in 2D and 3D porous media spanning a broad range of geometries.

We consider incompressible, inertia-free, viscoelastic flow through a porous medium [2], [?, ?]. The rate of work done by the imposed pressure drop, per unit volume of the medium, $\frac{\Delta P}{L}\cdot\frac{Q}{A}$, is balanced by the time ($t$)-and volume ($V$)-averaged rate of viscous dissipation, $\langle\boldsymbol{\tau}:\nabla\mathbf{u}\rangle_{t,V}$; here, $\boldsymbol{\tau}$ and $\mathbf{u}$ are the fluid stress tensor and velocity, respectively. Previous work [2] applied a Reynolds decomposition that separates this dissipation into mean and fluctuating components, $(\boldsymbol{\tau}:\nabla\mathbf{u})$ and $\langle\boldsymbol{\chi}\rangle_{t,V}\sim\langle\mathbf{s}^{\prime}:\mathbf{s}^{\prime}\rangle_{t,V}$, respectively, where $\mathbf{s}$ is the rate-of-strain tensor and a prime denotes the fluctuating part. The fluctuating component, $\langle\boldsymbol{\chi}\rangle_{t,V}$, captures the added viscous dissipation arising from unsteady flow fluctuations generated by the elastic instability. That decomposition, however, used a generalized-Newtonian (power-law) constitutive model of the fluid that neglects flow history-dependent polymeric stresses, and therefore omits the extensional contribution that we expect is required to fully capture flow thickening.

Instead, we extend this decomposition by separating the mean dissipation into its shear $(\tau_{ij,i\ne j})$ and extensional $(\tau_{ij,i=j})$ components. Rewriting the extensional stresses in terms of the extensional viscosity $\eta_E=\frac{\tau_{xx}-\tau_{yy}}{\dot{\epsilon}}$ and the extensional strain rate $\dot{\epsilon}=\frac{\partial u_y}{\partial x}$ yields an additional dissipation contribution $\langle\eta_E\dot{\epsilon}^{2}\rangle_{t,V}$, which captures the energy cost of stretching polymer chains under

The flow of viscoelastic polymer solutions through porous media at low Reynolds number ($\operatorname{Re} \ll1$) underlies diverse chemical production [1]–[4], additive manufacturing [5, 6], separation [7]–[10], energy production [11]–[18], and environmental remediation [19]–[22] processes. And yet, predicting how pore-scale dynamics set the macroscopic flow resistance, quantified by the “apparent viscosity” $\eta_{\mathrm{app}} \equiv\frac{\Delta P}{L}\frac{k}{Q/A}$, remains an open problem; here, Darcy’s law relates the total pressure drop $\Delta P$ across a length $L$ of the medium, the volumetric flow rate $Q$ through the cross section $A$, and the permeability $k$ [23]. Indeed, despite the typical shear-thinning rheology of these fluids, they anomalously ‘flow thicken’ in a porous medium: $\eta_{\mathrm{app}}$ increases, not decreases, above a threshold flow rate [24]–[32]. The underlying reason has been a mystery for over half a century. Various mechanisms have been proposed [2, 11, 13, 29]–[45], [47]–[?], [?]–[?], [?]–[?], [?]–[?], [?]–[?], but quantifying their relative contributions to flow thickening has been challenging due to the opacity of typical media. Instead, prior work has relied on numerical simulations to characterize pore-scale flow dynamics [51], [52], [?], [?], [?]–[?] or on semi-empirical models with adjustable fitting parameters [?]–[?]. flow. Combining all three contributions in Darcy’s law and normalizing by the bulk shear viscosity $\eta_{\mathrm{I}}$ evaluated at the characteristic pore-scale shear rate $\dot{\gamma}$ then gives our first main result:

where $\phi$ is the medium porosity and $D_p$ is the mean diameter of the pillars (2D) or grains (3D) composing the solid phase of the medium. The first term on the right-hand side of (1) represents the flow resistance from the bulk shear viscosity of the fluid. The second term reflects the viscous dissipation generated by unsteady flow fluctuations arising from an elastic instability when the Weissenberg number comparing elastic to viscous stresses in the fluid, $\mathrm{Wi}\equiv N_1(\dot{\gamma})/2\sigma(\dot{\gamma})$, exceeds a threshold value $\mathrm{Wi}_{\mathrm{c}}$ [2]; here, $N_1$ and $\sigma$ are the first normal stress difference and shear stress, respectively. The third term is the new extensional contribution, where the Trouton ratio $\mathrm{Tr}\equiv\eta_{\mathrm{E}}/\eta_0$ [8, 17] compares the Wi-dependent transient extensional viscosity to the zero-shear viscosity, $\eta_0$, and the characteristic extensional strain rate is $\dot{\epsilon}=\frac{Q\zeta}{\phi AD_p}$, where $Q/(\phi A)$ is the average flow speed through the pore space and $D_p/\zeta$ is the effective length scale over which that velocity changes as the fluid is squeezed past each obstacle; the value of the numerical factor $\zeta$ is determined from the unit-cell geometry [?].

To evaluate the accuracy of (1), we perform experiments using a polymer solution flowing through porous media of controlled ordered geometries Figure 1. The solution is composed of $c=300$ ppm ($\approx0.5c^*$, the critical overlap concentration) partially hydrolyzed polyacrylamide (18 MDa, 30% hydrolyzed) in a viscous solvent composed of 82.6% glycerol, 10.4% dimethyl sulfoxide, 6% ultrapure Millipore water, and 1% sodium chloride (by weight). This fluid is weakly shear-thinning and highly elastic, as quantified by shear rheology in Figure 1. We characterize the fluid’s extensional response using capillary break-up rheometry (CaBER); the Trouton ratio $\mathrm{Tr}$ grows with accumulated Hencky strain $\epsilon_{\mathrm{Hencky}}$ from the Newtonian value $\mathrm{Tr}_0=3$, saturating at $\mathrm{Tr}_{\infty}=1500$ as the stretching polymer chains approach their finite extensibility Figure 1 [8].

We study flow through two classes of ordered porous medium: (i) 2D millifluidic arrays of hexagonally arranged pillars, fabricated by stereolithographic 3D printing, in two configurations (‘staggered’ and ‘aligned’ relative to the mean flow direction), and (ii) 3D consolidated sphere packings, fabricated by selective laser-induced etching of fused silica, in two configurations (simple cubic, SC, and body-centered cuboid, BC) Figure 1. The 3D packings are rendered transparent by refractive-index matching of the solvent to the fused-silica matrix [?, ?], enabling direct flow visualization at the pore scale by confocal microscopy. We seed the polymer solution with 1 $\mu$m fluorescent tracer particles and measure the time-resolved velocity field $\mathbf{u}$ by particle-image velocimetry [16] in each unit cell across a range of imposed Weissenberg numbers. Simultaneously, we use a pressure transducer connected in parallel to measure the pressure drop across the medium [?], which we convert to $\eta_{\mathrm{app}}$ via Darcy’s law.

We first examine flow through the 2D hexagonal pillar arrays Figure 2. Bulk pressure drop measurements show that flow thickening sets in above $\mathrm{Wi}_{\mathrm{c}}\approx0.3$ in both configurations [points in Figure 2], coincident with the onset of the elastic instability [?]. However, the elastic instability contribution to the flow resistance [second term in (1)] falls far short of the measured resistance, as shown by the black curves in Figure 2, indicating that the extensional contribution is required. We evaluate this new contribution [third term in (1)] by extracting, from the measured pore-scale velocity fields, the probability density $p(\epsilon_{\mathrm{Hencky}})$ of accumulated Hencky strain experienced by fluid elements along Lagrangian trajectories integrated over one polymer relaxation time, mapping it to a distribution $p(\mathrm{Tr})$ of local Trouton ratios via the CaBER-measured $\operatorname{Tr}(\epsilon_{\mathrm{Hencky}})$ [][(#fig:1b)], and computing the local flow-field-weighted average $\langle\operatorname{Tr}\rangle_{t,V}=\int_{\operatorname{Tr}_0}^{\operatorname{Tr}_\infty}\operatorname{Tr}\,p(\operatorname{Tr})\,\mathrm{d}\operatorname{Tr}$ [?].

With the extensional contribution included, (1)—evaluated using only fluid rheology and pore-scale flow field measurements—quantitatively recovers the measured flow thickening in both 2D arrays, as shown by the solid colored curves in fig:2c,d [?].

To test the generality of (1), we repeat this analysis for flow through the 3D consolidated sphere packings [fig:3a,b]. Flow thickening again sets in coincident with the elastic instability, here at $\mathrm{Wi}_c\approx2-3$ [points in fig:3c,d], but the elastic instability term is again insufficient to account for the measured resistance [black curves]. Instead, just as in the 2D case, adding the extensional contribution—evaluated from the measured pore-scale velocity fields—recovers the measured resistance in both packings [solid colored curves] [?]. Altogether, the agreement of (1) with our experiments across 2D and 3D geometries, spanning porosities $\phi$ from 0.28 to 0.85, establishes that polymer extension is the dominant additional contribution required to capture flow thickening

in ordered porous media beyond the shear and elastic-instability terms.

Where in the pore space does this extension manifest?

One might expect pore constrictions (so-called “throats”) to be the locations of strongest polymer stretching: each throat compresses streamlines into a smaller cross section, generating an extensional velocity gradient that stretches polymer chains as they pass through [39, ?, ?]. However, the maximum accumulated Hencky strain that throats alone can generate in each unit cell is $\epsilon_{\mathrm{Hencky}}\approx0.3-1.5$ [?]—much smaller than the measured strains, which routinely exceed the polymer coil-stretch threshold of $\epsilon_{\mathrm{Hencky}}\approx2-3$ [28, 29]. Instead, prior work [11, 40] demonstrated that stagnation points in the pore space, not throats, generate the strongest local velocity gradients and therefore most effectively stretch polymer chains [?, ?, ?, ?]. If extension is stagnation point-dominated, then each stagnation point should contribute to flow thickening additively, and the macroscopic extensional resistance should scale linearly with the number of stagnation points per unit cell, $n_{\mathrm{SP}}$.

We test this hypothesis by measuring the resistance from polymer extension at a single stagnation point and then multiplying by $n_{\mathrm{SP}}$, counted a priori from the geometry of each medium [?]. We use an optimized cross-slot rheometer [41, 42] as a single-stagnation point “model pore”: in extensional flow mode, one free stagnation point sits at the device center, generating a planar extensional flow whose excess pressure drop relative to the corresponding shear-only configuration, $\Delta P_{\mathrm{ex}}/\Delta P_{\mathrm{sh}}$, isolates the per-stagnation point contribution of polymer extension [?]. Decomposing the Darcy pressure drop into shear and extensional parts then yields our second main result:

As shown by the dashed colored curves in Figs. 2(c,d) and 3(c,d), Eq. (2) again recovers the measured flow resistance across all the 2D and 3D geometries tested, with all input parameters measured independently. This agreement between the theory and experiment establishes that the macroscopic flow resistance of viscoelastic polymer solutions in ordered porous media is set, additively, by polymer extension at stagnation points.

Our experiments so far have focused on ordered porous media. We next test whether the framework extends to disordered porous media as well using experiments in two model 3D geometries: lightly sintered packings of near-monodisperse glass beads or of crushed glass, both designed to be refractive-index matched, enabling in situ flow visualization. As detailed in [?], in both cases, the elastic instability contribution [second term of (1)] captures the resistance near the onset of flow thickening but progressively underpredicts the measured increase in $\eta_{\mathrm{app}}$ at higher Wi—suggesting that polymer extension also contributes to flow thickening in disordered media. Our analysis confirms this idea: (1), evaluated using the experimentally-measured Hencky-strain distribution, captures the macroscopic resistance across the full Wi range in the disordered packings Figure 4—again with $D_{\mathrm{p}}/\zeta\approx0.6-0.8D_{\mathrm{p}}$ characteristic of an extensional length scale set by the local pore geometry. Our framework therefore captures flow thickening across porous media of different geometries.

We have focused on shear, elastic instability-generated fluctuation, and extensional contributions to flow thickening in porous media, and in doing so have necessarily made three key simplifications. First, we assumed that shear-generated normal-stress differences contribute negligibly, as supported by an upper-bound estimate combining the bulk measurements of $N_1$ with the measured pore-scale velocity-gradient field [?]; incorporating this contribution explicitly [53] is a natural extension of our framework. Second, our continuum description of the fluid neglects polymer-wall interactions, justified here because the polymer size (radius of gyration $R_{\mathrm{g}} \approx200$ nm) is more than three orders of magnitude smaller than the pore throats ($\gtrsim300$ $\mu$m). These interactions will likely become more important in more confined media [38, 41]. Finally, we have not included the additional contribution $\boldsymbol{\tau}_{\mathrm{p}}:\nabla\mathbf{u}$, where $\boldsymbol{\tau}_{\mathrm{p}}$ is the polymer stress tensor, that accounts for the storage and dissipation of elastic energy during polymer deformation and relaxation [21, 26, ?], which could be quantified using birefringence [41, ?, ?, ?, ?] or by reconstructing polymer stress fields along Lagrangian fluid trajectories [23, 24, 25], both of which remain active research frontiers. Even with these simplifications, our framework provides a quantitative pore-to-macroscale link for flow thickening across both ordered and disordered porous media—which could thus guide efforts to engineer targeted flow behavior via fluid rheology [?, ?] and pore geometry [?, ?].

Acknowledgments

We thank H. A. Stone and R. J. Poole for insightful discussions. The authors acknowledge the use of the Imaging and Analysis Center (IAC) operated by the Princeton Materials Institute at Princeton University, which is supported in part by the Princeton Center for Complex Materials (PCCM), a National Science Foundation (NSF) Materials Research Science and Engineering Center (MR-SEC; DMR-2011750). S.S.D. also acknowledges support from the Camille Dreyfus Teacher-Scholar Program of the Camille and Henry Dreyfus Foundation. S.J.H. and A.Q.S. also acknowledge financial support from the Japanese Society for the Promotion of Science (JSPS, Grant Nos. 24K07332, and 24K00810).

Polymer solution

To prepare the polymer solution, we gently mix a stock solution of 5000 ppm (by weight) partially hydrolyzed polyacrylamide (HPAM, 18 MDa, 30% hydrolysis; Polysciences) in ultrapure Millipore water on a rotor mixer for 24 hours. We then dilute this stock solution to 300 ppm HPAM with the remaining solvent components (by weight): 82.6% glycerol (ThermoFisher Scientific), 10.4% dimethyl sulfoxide (Sigma-Aldrich), and 1% sodium chloride (Sigma-Aldrich) and mix for another 24 hours on a stir plate at 60 rpm to avoid mechanical degradation of the polymer. The solution is then allowed to rest for 24 hours before any rheology measurements or flow experiments. All polymer solutions are used within 1 month of preparation.

The addition of excess salt (170 mM NaCl) sufficiently screens the polyelectrolyte charges such that HPAM is expected to adopt a flexible coil configuration in solution [1]. For this solvent, the overlap concentration at which polymer intermolecular interactions occur is $c^{*} = 600$ ppm as determined from shear rheology of a dilution series to determine the intrinsic viscosity [2]; thus, the polymer solution is in the dilute regime, $c/c^{*} = 0.5$. From the overlap concentration $c^{*} = \frac{M_{W}}{\frac{4}{3}\pi R_{g}^{3}N_{A}}$, we estimate the radius of gyration as $R_{g} = 220$ nm, where $M_{W}$ is molecular weight and $N_{A}$ is Avogadro’s number. The fully extended length of the polymer chain is given by the contour length, $L_{c} = NL_{\mathrm{mon}} \approx40$ μm with $N$ the degree of polymerization and $L_{\mathrm{mon}}$ the monomer length.

Fluid rheology

We measure the steady shear rheology of the polymer solution using a stress-controlled Anton Paar MCR501 rheometer with a truncated cone-plate geometry (CP50-2: 50 mm diameter, $2^\circ$, 53 μm gap), temperature-controlled at $25^\circ$C. We measure steady-shear flow curves of the shear stress, $\sigma(\dot{\gamma})$, and first normal stress difference, $N_{1}(\dot{\gamma})$, over shear rates of $\dot{\gamma} = 0.1-30$ s$^{-1}$. For shear rates $\dot{\gamma} > 30$ s$^{-1}$, we observe an elastic instability in the cone-plate geometry marked by an apparent thickening of the shear stress [3, 4, 5, 6]. The minimum resolvable $N_{1}$ is $\approx10$ Pa due to the normal force sensitivity of the rheometer. We perform flow curve measurements for three replicates with fresh samples for each sweep and report the standard deviation as error bars. We fit a power law model to the shear stress, $\sigma(\dot{\gamma}) = K\dot{\gamma}^{n}$, and first normal stress difference, $N_{1}(\dot{\gamma}) = K_{N_{1}}\dot{\gamma}^{n_{N_{1}}}$, where $K$ and $K_{N_{1}}$ are material prefactors and $n$ and $n_{N_{1}}$ are power law indices. For the dilute HPAM solution, $K = 0.323 \pm0.001$ Pa·s$^{n}$, $n = 0.93 \pm0.1$, $K_{N_{1}} = 0.81 \pm0.05$ Pa·s$^{n_{N_{1}}}$, and $n_{N_{1}} = 1.43 \pm0.03$. In the absence of a low shear rate viscosity plateau for power law fluids, we approximate the zero-shear viscosity $\eta_{0}$ as $\eta|_{\dot{\gamma}=0.1\ \mathrm{s}^{-1}} \approx340$ mPa·s.

We characterize the transient extensional flow behavior of the polymer solution using a Haake capillary breakup extensional rheometer (CaBER). In this method, a small volume of fluid is loaded between two plates at an initial plate separation of 1 mm, the plates are separated at 1 m/s to a final separation of 6 mm to create a thin filament, and the resulting thinning of the filament diameter is recorded over time [7, 8]. In the elasto-capillary regime, filament thinning is governed by a balance of elastic stresses resisting break-up and the capillary forces promoting pinch-off [9]. The Hencky strain over the course of the measurement is $\epsilon_{\mathrm{Hencky}} = 2\cdot\ln(D_{0}/D_{\mathrm{mid}}(t))$, where $D_{0}$ is the plate diameter and $D_{\mathrm{mid}}(t)$ is the filament diameter at the mid-plane. An analysis procedure for balance on the filament [7, 9] enables determination of the Trouton ratio (Tr; ratio of transient extensional viscosity to zero-shear viscosity $\eta_{E}/\eta_{0}$) given that the surface tension and density of the fluid are known. For our polymer solution, we measure a surface tension of 61.2 mN/m (Biolin Attension Theta; pendant drop method) and density of 1234 kg/m$^{3}$ (Anton Paar DMA 35 densitometer). We measure five replicates in the CaBER device at $25^\circ$C and find good reproducibility. We fit the measured Trouton ratio versus Hencky strain data with a logistic function to capture the low-strain Newtonian-like behavior and high-strain finite extensibility of the polymer:

where $\mathrm{Tr}_0 = 3$ is the low-strain Newtonian limit, $\mathrm{Tr}_{\infty} = 1500$ is the finite extensibility limit, $k$ is the growth rate, and $\epsilon_{H,0.5}$ is the Hencky strain at which the Trouton ratio reaches its half-maximum. Here, we determine $\mathrm{Tr}_{\infty} = 1500$ based on the upper limit of the elastocapillary regime. We use the resulting functional form to convert Hencky strains measured from the flow field to estimates of the “local” apparent Trouton ratio (i.e. a pore-scale estimate of transient extensional viscosity).

Porous medium fabrication

2D hexagonal pillar arrays

We fabricate millifluidic 2D arrays consisting of pillars of diameter $D_p = 600$ $\mu$m arranged in a hexagonal unit cell with center-to-center spacing of $l_p = 1.44$ mm. These dimensions correspond to a ratio of $l_p/D_p = 2.4$, which we have chosen to reproduce following the devices used by Walkama et al. [10] and Haward et al. [11]. In the staggered configuration, the streamwise pillar-pillar spacing is $l_s = 2.494$ mm, while for the aligned array, $l_s = l_p$. The difference in $l_s$ for the two geometries is a result of the imposed flow direction relative to the unit cell orientation. The medium porosities are $\phi_{\mathrm{stag}} = 0.848$ and $\phi_{\mathrm{aligned}} = 0.836$. We determine the Darcy permeability of each medium, $k$, by measuring the pressure drop across the device for the flow of water applied in a stepped flow rate ramp sequence and using Darcy’s law. The resulting Darcy permeabilities are $k_{\mathrm{stag}} = 895 \pm20$ $\mu$m$^2$ and $k_{\mathrm{aligned}} = 929 \pm19$ $\mu$m$^2$. The arrays have lengths of $L_{\mathrm{stag}} = 33$ mm and $L_{\mathrm{aligned}} = 32.3$ mm and widths of $W = 8.4$ mm. The channel depth is 1 mm. In total, the staggered array contains 149 pillars and the aligned array contains 157 pillars.

We design the channel using CAD software (Onshape) and 3D-print the main body using a FormLabs Form 3 stereolithography 3D printer with a proprietary clear resin consisting of methacrylate oligomers and photoinitiators (FLGPCL04). We assemble the full millifluidic device by sandwiching a $\sim1$ mm layer of polydimethylsiloxane (PDMS; Dow SYLGARD 184, 8.5:1.5 by weight base to curing agent ratio) between the 3D printed base and a laser-cut transparent acrylic top sheet (McMaster-Carr) fitted with screwholes. The PDMS layer ensures a watertight seal while allowing for optical access within the channel. The device is assembled by tightly screwing the components together. We glue flexible Tygon tubing (McMaster-Carr) into the inlet and outlet of each device using a watertight 2-part epoxy (JB Marine-Weld).

3D ordered sphere packings

The fabrication protocol for the 3D consolidated sphere packings is described in previous work [12]. The bead diameter in both packings is $D_p = 780$ $\mu$m, and the diameter of the cylindrical consolidations is $D_{\mathrm{neck}} = 170$ $\mu$m ($D_{\mathrm{neck}}/D_p = 0.22$). Neglecting the volume of the consolidations, the simple cubic (SC) unit cell has dimensions of $D_p \times D_p \times D_p$. Similarly, the body-centered cuboid unit cell has dimensions of $D_p \times D_p \times\sqrt{2}D_p$. Note that the BC geometry is a rectangular cuboid unit cell rather than an idealized cubic body-centered unit cell. The medium porosities are $\phi_{\mathrm{SC}} = 0.476$ and $\phi_{\mathrm{BC}} = 0.284$. The medium Darcy permeabilities are $k_{\mathrm{SC}} = 567 \pm9$ $\mu$m$^2$ and $k_{\mathrm{BC}} = 197 \pm2$ $\mu$m$^2$. The total lengths of each packing are $L_{\mathrm{SC}} = 3.84$ cm and $L_{\mathrm{BC}} = 2.74$ cm, and the cross-sectional area is $A = 0.16$ cm$^2$.

Flow parameters

For the 2D pillar arrays, we define the characteristic shear rate as $\dot{\gamma} = U/R_p$, where $U = Q/(HW\phi)$ is the average flow velocity and $R_p$ is the pillar radius. For the 3D sphere packings, we define the characteristic interstitial shear rate as $\dot{\gamma}_I = Q/(A\sqrt{k\phi})$ [13, 14, 15] and characteristic interstitial velocity as $U_I = Q/(A\phi)$.

We define the Reynolds number, the ratio of inertial to viscous stresses, as $\mathrm{Re} = \rho U D_p/\eta_I$, where $\rho$ is the fluid mass density, $U$ is the characteristic flow speed, and $\eta_I$ is the bulk shear viscosity evaluated at the characteristic shear rate. In our experiments, $\mathrm{Re} < 10^{-2}$ which renders the effects of inertia negligible. Thus, the onset of flow instability is assumed to arise solely from fluid elasticity. Additionally, given the weak shear-thinning nature of the polymer solution ($n = 0.93$), we do not expect shear-thinning to influence the flow in our experiments.

The Weissenberg number quantifies the ratio of elastic to viscous stresses in the fluid, which we define here as $\mathrm{Wi} = N_1(\dot{\gamma})/[2\sigma(\dot{\gamma})]$. For a given imposed flow rate $Q$, we determine a nominal imposed $\mathrm{Wi}_I$ evaluated using the fluid shear rheology and characteristic shear rate.

Flow imaging and pressure drop measurements

We prepare each experiment by first flushing the pristine media with ultrapure Millipore water to saturate the pore space and remove any air bubbles. In the 2D devices, we sequentially flush with Millipore water and pure glycerol before pre-filling the medium with the polymer solution at a constant low flow rate of $Q = 1$ mL/hr for at least 3 hours using a Harvard Apparatus PHD2000 syringe pump. In the 3D devices, we sequentially flush with Millipore water, isopropanol, Millipore water, and pure glycerol before pre-filling the medium with the polymer solution at a constant low flow rate of $Q = 1$ mL/hr for at least 3 hours. This procedure ensures that the pore space is entirely saturated with the polymer solution prior to the flow experiment.

Flow imaging. For each flow imaging experiment, we apply a stepped flow rate ramp spanning a range of $\mathrm{Wi}$ to probe steady and unsteady flow conditions. After each step change in flow rate, we let the flow equilibrate for at least one hour before recording any data. To visualize the flow, we seed the fluid with 1 µm carboxylate-modified, fluorescent, polystyrene tracer particles (Invitrogen; yellow-green) that have excitation at 480-510 nm and emission at 515 nm. The tracer particle seeding concentrations for the 2D and 3D flow experiments are 10 ppm and 7 ppm, respectively. We place the medium on the stage of a Nikon A1R+ laser scanning confocal fluorescence microscope, where particles are excited by a 488 nm laser and detected with a 500-550 nm sensor. We record 2 min videos per imaging location at 30 fps using the confocal resonant scanner with a raster scanning path.

In the 2D arrays, we use a 4× objective to interrogate a 3166.85 µm x 3166.85 µm field of view with optical section thickness 36.95 µm and pixel size 6.19 µm/px (512 x 512 image). We focus the imaging plane at the channel half-height to avoid wall effects.

In the 3D packings, we use a 10× objective to interrogate a 1272.79 µm x 1272.79 µm field of view with optical section thickness 7.39 µm and pixel size 2.46 µm/px (512 x 512 image). To capture the $z$-variations in the flow field, given the 2D flow imaging constraints of confocal microscopy, we image nine selected $z$-planes spanning the full unit cell ($\delta z = 50$ µm) corresponding to depths ranging from ∼ 400 − 1200 µm from the bottom of the medium. We aggregate the measured flow fields across $z$ at each $\mathrm{Wi}$ for analysis.

In each device, we image multiple unit cells located at various positions throughout the pore space in each geometry. Sufficiently far from the inlet (∼ 1-2 unit cells) and walls, we do not observe spatial variations in the flow field at lengthscales larger than the characteristic unit cell. Thus, for all subsequent analysis, we analyze flow fields corresponding to one unit cell located in the middle of each porous medium.

Image processing. We use a MATLAB particle image velocimetry toolbox (PIVlab) [16] to obtain time-resolved, 2D velocity fields with a temporal resolution of 33 ms. In PIVlab, we use the FFT window deformation algorithm with 4 passes of decreasing interrogation area (64-32-16-16 pixels) and 50% overlap. The spatial PIV grid resolutions are 50.3 µm (2D arrays) and 20.2 µm (3D packings).

We generate pathline images and videos of the flow using ImageJ by averaging the mean pixel intensity over successive frames.

Pressure drop measurements. We quantify the macroscopic flow behavior using pressure drop measurements obtained by connecting a pressure transducer (Omega; 0-30 psi) in parallel to the inlet and outlet of each device. We impose a stepped sequence of increasing flow rate using a syringe pump, allowing the flow at each flow rate to equilibrate for at least one hour. We then time-average the pressure drop signal (10 Hz sampling rate) after the equilibration period and adjust for measurement noise by subtracting the pressure drop for no-flow conditions. Figure S1 presents the raw pressure drop-flow rate measurements for the viscoelastic polymer solution in the four porous media. The solid line corresponds to Darcy’s law evaluated using the zero-shear viscosity, $\eta_0$. In all devices, the pressure drop follows the linear Darcy prediction at low flow rates; however, above a threshold flow rate, a pronounced positive deviation from Darcy’s law is observed—a “flow-thickening” of the polymer solution. We convert the pressure drop measurements into a dimensionless measure of flow resistance by calculating an apparent viscosity from Darcy’s law and normalizing by the bulk shear viscosity measured from rheology, $\frac{\eta_{\mathrm{app}}}{\eta_l} = \frac{k\Delta P}{L(Q/A)\eta_l}$. A value of $\frac{\eta_{\mathrm{app}}}{\eta_l} = 1$ indicates that the apparent viscosity of the fluid in the porous media is equivalent to the bulk shear viscosity as measured in a cone-plate shear rheometer; values of $\frac{\eta_{\mathrm{app}}}{\eta_l} > 1$ correspond to a higher viscosity than predicted by bulk shear rheology alone. We do not observe any enhanced porous medium drag reduction effects ($\frac{\eta_{\mathrm{app}}}{\eta_l} < 1$) for this HPAM solution in our experiments.

Derivation of power balance incorporating an apparent extensional viscosity

To derive the dimensionless flow resistance model incorporating the local apparent extensional viscosity, we start from the mechanical energy balance obtained by dotting the Cauchy momentum equation with velocity:

where $\mathbf{u}$ is velocity, $P$ is pressure, $\boldsymbol{\tau}$ is the fluid stress tensor, and $\mathbf{g}$ is gravitational acceleration. Following the derivation by Browne and Datta [2], we neglect the inertial kinetic energy terms ($\operatorname{Re} \ll1$), assume incompressible flow ($\nabla\cdot\mathbf{u} = 0$), neglect viscous work at the walls of the porous medium, and neglect gravitational work. Integrating the resulting energy balance and averaging over volume and time, we obtain a macroscopic energy balance,

which balances the pressure work done by the fluid to the energy dissipation in the fluid. Employing Darcy’s law, this expression can be recast in terms of the apparent viscosity in a porous medium:

where we have normalized by the bulk shear viscosity of the fluid to quantify a dimensionless measure of the flow resistance.

Previously, Browne and Datta [2] used a Reynolds decomposition for the viscous dissipation term, $\langle\boldsymbol{\tau} : \nabla\mathbf{u} \rangle_{t,V}$, to separate the mean and fluctuating flow components. In doing so, they derive the following power balance that accounts for resistances from the shear viscosity and unstable viscous dissipation arising from an elastic instability:

Here, $\langle\chi\rangle_{t,V} \sim\langle\mathbf{s}' : \mathbf{s}' \rangle_{t,V}$ is the unstable viscous dissipation directly measured from the pore-scale flow fields and $\mathbf{s}'$ is the fluctuating components of the rate-of-strain tensor. The growth in the unstable viscous dissipation is described by a power law: $\langle\chi\rangle_{t,V} = \alpha(\mathrm{Wi}-\mathrm{Wi}_{c})^{\beta}$, where $\alpha$ and $\beta$ are fitting parameters and $\mathrm{Wi}_{c}$ is the critical Wi for the onset of instability. This derivation assumes a generalized Newtonian fluid model—specifically, a power law fluid—which captures the non-Newtonian shear-thinning behavior of polymer solutions, but neglects any elasticity or enhanced extensional viscosity beyond that of a Newtonian fluid ($\mathrm{Tr}_{\mathrm{Newt}} = 3$ for uniaxial extension [17, 18, 19]). A large increase in extensional viscosity can be observed for elastic polymer solutions resulting from the alignment and stretching of flexible polymer chains under flow [20, 21, 22]. The increase in extensional viscosity directly results from an increase in the polymeric tensile stresses, $\tau_{p,ii}$, which are strongly flow history-dependent. The relaxation of polymer stresses over finite timescales imparts the fluid with a “memory” of its deformation history. Experimentally, it is challenging to directly resolve these tensile stresses $\tau_{p,ii}$ and this problem remains a subject of ongoing research [23, 24, 25]. Thus, we adopt an alternative approach to evaluating the enhanced extensional stresses through a local apparent extensional viscosity rather than attempting to evaluate the polymer stresses directly.

Returning to the decomposition of the viscous dissipation term proposed by Browne and Datta [2]:

we now further decompose the mean dissipation term to separate the shear and extensional stresses, producing,

The terms in the first bracket on the right side of (5) contain the tensile stresses and the terms in second bracket contain the shear stresses. Here, we assume that the terms containing shear stresses are adequately captured by the shear viscosity, $\eta_I(\dot{\gamma})$, from bulk rheology. The extensional stress terms can be rewritten as,

The term $(\tau_{yy}-\tau_{zz})$ is the second normal stress difference, $N_2$, which is typically negligible for dilute polymer solutions [26, 27]; hence, we neglect the second term on the right side of (6). Finally, we introduce the extensional viscosity, $\eta_E = \frac{\tau_{xx}-\tau_{yy}}{\dot{\epsilon}}$, and strain rate, $\dot{\epsilon} = \frac{\partial u}{\partial x}$ to arrive at the modified power balance:

The last term in (7) now accounts for enhanced extensional stresses arising from polymer stretching quantified via an apparent extensional viscosity.

Hencky strain and determination of local Trouton ratio

Measurement of Hencky strain

Our approach to evaluating the extensional resistance term in (7) first requires measuring distributions of pore-scale accumulated Hencky strains which quantifies the net extensional deformation of polymers. We determine the Hencky strain for a fluid particle by initializing the particle at a given location within a pore and integrating forward in time using the experimental velocity field over a total interval equal to the longest relaxation time of the polymer; the longest relaxation time of the 300 ppm HPAM solution is $\lambda_0 = 0.5$ s from capillary break-up rheometry. The longest relaxation time corresponds to the maximum time to which polymeric elastic stresses can persist before being dissipated by relaxation mechanisms. Mathematically, the Hencky strain is,

which we numerically approximate as,

Here, $\lvert\dot{\epsilon} \rvert= \left\lvert\frac{\partial u}{\partial x} \right\rvert+ \left\lvert\frac{\partial v}{\partial y} \right\rvert$ is the local strain rate and $\delta t$ is the time step set by the temporal resolution of the experimental flow field. We use the absolute value of the strain rate to track the cumulative extensional deformation over stretching events. For each time step, we evaluate $\lvert\dot{\epsilon} \rvert$ at the current fluid particle position and use the corresponding fluid velocity $\mathbf{u}(\mathbf{x},t)$ to advance the particle’s position to the next time step, $\delta\mathbf{x} = \mathbf{u}(\mathbf{x},t)\delta t$. In total, we repeat this Hencky strain calculation for 15 randomly selected time points in the imaging window and for 5 randomly chosen initial particle positions to generate a probability density function (p.d.f) of accumulated Hencky strains for a given imposed Wi$_{\mathrm{I}}$. Note that the limited imaging field-of-view sets a spatial constraint on the maximum distance a fluid particle may travel over the relaxation time integration interval.

Figure S2(a) presents the inverse cumulative probability density functions (1-CDF) of Hencky strains in the four porous medium geometries. Generally, the magnitude of accumulated Hencky strain increases with increasing Wi$_{\mathrm{I}}$. Importantly, we observe that the Hencky strains can exceed values of $\epsilon_{\mathrm{Hencky}} = 2$–$3$, which is a benchmark threshold for the polymer coil-stretch transition [28, 29]. Above the coil-stretch transition, velocity gradients are sufficiently large enough to orient and stretch polymer chains from their equilibrium coil configurations, producing a concomitant increase in the fluid extensional viscosity relative to the shear viscosity [22, 30].

Weighted Trouton ratio calculation

Using the distributions of measured Hencky strains, we can further determine distributions of local apparent Trouton ratios by using the relation $\operatorname{Tr}(\epsilon_{\mathrm{Hencky}})$ in (3) from capillary break-up rheometry. Note that the specific parameters in the logistic equation for $\operatorname{Tr}(\epsilon_{\mathrm{Hencky}})$ are fluid-specific. For each imposed $\mathrm{Wi}_{\mathrm{I}}$, we obtain a probability density distribution of local apparent Trouton ratios, $p(\operatorname{Tr})$, describing fluid particles that are advected and deformed by the flow.

Finally, to evaluate the extensional resistance term in (7), we first substitute for the transient extensional viscosity, $\eta_{\mathrm{E}}=\eta_{0}\operatorname{Tr}$, and nominal strain rate, $\dot{\epsilon}=\frac{Q\zeta}{A\phi D_{\mathrm{p}}^{2}}$, where $\zeta$ is a dimensionless, geometry-dependent parameter (see §). The resistance term can be rewritten as,

We use the probability density distributions, $p(\operatorname{Tr})$, to calculate $\langle\operatorname{Tr}\rangle_{t,V}$ by calculating the probability density-weighted expectation value using numerical integration:

We determine $\langle\operatorname{Tr}\rangle_{t,V}$ for each imposed $\mathrm{Wi}_{\mathrm{I}}$ in the flow experiments and linearly interpolate or extrapolate between these discrete points to evaluate the power balance model; Figure S2 shows the flow field-weighted Trouton ratios as a function of $\mathrm{Wi}_{\mathrm{I}}$ in the four porous medium geometries. Open circles are the experimental values and the solid curves are linear interpolation/extrapolations. Physically, $\langle\operatorname{Tr}\rangle_{t,V}$ can be interpreted as a medium-averaged value of the local apparent Trouton ratio weighted by the pore-scale flow field—this weighting captures variability in the extensional flow behavior within the pore space as different fluid particles experience different flow kinematics depending on their Lagrangian trajectory.

Theoretical Hencky strain due to contractions set by geometry

A useful feature of the ordered porous medium geometries in this work is that the unit cell void spaces can be described mathematically using geometrical considerations. As such, we can analytically derive the predicted Hencky strains for extension caused by contractions in the pore space alone (i.e. neglecting extensional flow at stagnation points). Given the position-dependent cross-sectional area of the void space, $A(x)$, the position-dependent Hencky strain is, $\epsilon_{\mathrm{Hencky}}(x)=\ln\left(\frac{A_{\max}}{A(x)}\right)$. Here, $x$ denotes the streamwise position along the length of the unit cell. For the 2D arrays, the Hencky strain is equivalent to $\epsilon_{\mathrm{Hencky}}=\ln\left(\frac{L_{\max}}{L(x)}\right)$, where $L(x)$ is the length of the void space at a given position $x$ obtained by trivially dividing $A(x)$ by the constant channel height. Equations (S7)–(S10) below present the mathematical expressions for the void space in each unit cell. Here, $l_{\mathrm{p}}$ is the pillar center-center distance, $l$ is the streamwise length of the BC unit cell ($l=\sqrt{2}D_{\mathrm{p}}$) Figure S4, and $D_{\mathrm{p}}$ is the pillar (2D) or grain diameter (3D).

Figures S3(a, b) show schematics of the 2D pillar array unit cell geometries and characteristic length scales. Note that here we use a truncated rectangular unit cell rather than the full hexagonal unit cell for ease of calculation. Figures S3(c, d) present the local Hencky strain profiles for the staggered and aligned unit cells. Yellow stars mark the maximum values of Hencky strain: $\varepsilon_{\mathrm{Hencky}}^{\max}(\mathrm{stag.})=0.54$ and $\varepsilon_{\mathrm{Hencky}}^{\max}(\mathrm{aligned})=0.28$. This is a key result: the measured Hencky strains from the experimental flow fields Figure S2 are much larger than the maximum Hencky strains associated with contraction-generated extensional flow [Figure S3](c, d); thus, we hypothesize that extensional flow in the 2D arrays is dominated by stagnation points rather than contractions in the pore space geometry—ultimately motivating the resistance-per-stagnation point model in §.

Additionally, we use the analytical expressions of the void space cross-sectional area to estimate a position-dependent strain rate in each unit cell following Haward and Odell [31, 32]. For constant imposed 2D volumetric flux $q$, conservation of mass requires $q=\bar{V}\bar{L}=V(x)L(x)$, where $V(x)$ is the position-dependent flow speed and an overbar denotes a spatially-averaged quantity. Rearranging, we find $V(x)=q/L(x)$ and differentiation produces the local strain rate, $\dot{\varepsilon}(x)=\frac{dV(x)}{dx}$. The resulting strain rate profiles (non-dimensionalized by $q/D_p^2$) are shown in [Figure S3](e, f). Because of the circular nature of the pillars in this geometry, there is a numerical artifact where the analytical expressions for $\dot{\varepsilon}(x)$ diverge, denoted by the dotted lines on the plots.

We perform a similar analysis for the 3D sphere packing unit cells [Figure S4](a, b), now using the full void space cross-sectional area, $A(x)$. Figures S4(c, d) present the local Hencky strain profiles with maximum values of $\varepsilon_{\mathrm{Hencky}}^{\max}(\mathrm{SC})=1.54$ and $\varepsilon_{\mathrm{Hencky}}^{\max}(\mathrm{BC})=0.49$. Again, these maximum values are smaller than the Hencky strain values we measure from the flow fields in Figure S2, and we conclude that extensional flow in the 3D packings is dominated by stagnation points rather than contractions in the pore space geometry. Figures S4(e, f) present the local strain rate profiles in each unit cell non-dimensionalized by $Q/D_p^3$, determined from $\dot{\varepsilon}=\frac{dV(x)}{dx}$ and $V(x)=Q/A(x)$ with $Q$ a constant volumetric imposed flow rate. Yellow stars mark the locations of maximum strain rate. In this analysis, we neglect the consolidations to simplify the geometric expressions for the void space; however, we do not expect our conclusions regarding the contraction-generated versus stagnation point-generated extensional strains to change substantially with the inclusion of the consolidations given their small contribution to the total solid volume.

Physical interpretation of effective strain rate

As mentioned in §, we define the nominal strain rate in the porous medium as,

where we use $\zeta$ in (8) as a fitting parameter in the power balance model (7). Similar approaches have been used to define an effective interstitial shear rate [13, 14, 15] and nominal extensional strain rates [33, 34, 35, 36] in prior literature, where parameters analogous to $\zeta$ are used as geometry-dependent fitting parameters. However, previous analyses, to our knowledge, do not provide a physical interpretation of the corresponding fitting parameters. Here, we provide a physical interpretation and rationalization of the identified $\zeta$ values in the extensional viscosity-modified power balance model for the four ordered porous medium geometries.

Table S1 lists the experimental $\zeta$ values from the power balance models in each of the four porous medium geometries. A nominal strain rate can be approximated as a streamwise velocity difference over a characteristic length scale of extension: $\frac{\partial u}{\partial x}\approx\frac{\Delta u}{\Delta l}$. Our definition implies a velocity difference of $\Delta u\approx(U_I-0)=\frac{Q}{\phi A}$ over a length scale $\Delta l\approx D_p/\zeta$, where $U_I$ is the characteristic interstitial velocity. Table S1 summarizes the resulting extensional length scales $\Delta l\approx D_p/\zeta$ proposed by our analysis.

Geometry	Theory	Expt.	$1/\zeta$	$\alpha$
Stag. (2D)	0.443	0.53	1.89	0.72
Aligned (2D)	0.275	0.47	2.13	0.18
SC (3D)	0-5.52	1.13	0.89	—
BC (3D)	0-4.25	1.55	0.65	—

In the 2D arrays, the extensional length scale in both configurations is larger than the pillar diameter $D_p$ but smaller than the streamwise pillar spacing $l_s$. Previous work studying viscoelastic flows in cylindrical arrays has suggested that the relevant length scale of extension depends on the density of pillars (or solidity): small diameter pillars ($D_p/l_s\ll1$) act as “point sources” of stagnation points producing a dominant extensional length scale equal to the distance between point sources $l_s$, while extension in arrays of larger pillars with small spacing ($D_p/l_s\lesssim1$) is controlled by $D_p$ [36, 37, 38, 39]. Our results align with these prior studies, suggesting that the appropriate extensional length scale for the hexagonal arrays used in this work fall somewhere between these two extremes. Indeed, we may consider $D_p/\zeta$ as a linear combination of the two length scales $D_p$ and $l_s$ given by, $\frac{D_p}{\zeta} = \alpha D_p + (1-\alpha)l_s$, with $\alpha$ a weighting parameter ranging from $0-1$. We find $\alpha_{\mathrm{stag}} = 0.72$ and $\alpha_{\mathrm{aligned}} = 0.18$. A lower value of $\alpha$ (smaller weight on $D_p$ compared to $l_s$) in the aligned array compared to the staggered array can be rationalized by the increased screening of stagnation points with increasing Wi, as proposed by Haward et al. [11].

In the 3D packings, the values of $D_p/\zeta$ suggest that the relevant extensional length scale is smaller than the grain diameter $D_p$. Intuitively, this aligns with our previous work where we find that elastic instabilities in these consolidated sphere packings are generated from stagnation points arising at grain consolidations rather than stagnation points on the grain surfaces themselves [40]. As a result, we expect that the effective length scale of extension in the 3D ordered sphere packings lies between the diameter of the consolidations $D_{\mathrm{neck}}$ and the diameter of the grains $D_p$. For the specific geometries in this work, $D_{\mathrm{neck}}/D_p \approx0.22$, satisfying our hypothesis that $D_{\mathrm{neck}}/D_p \le1/\zeta\le1$.

Theoretical prediction of $\zeta$

The local strain rate analysis in § also enables theoretical prediction of $\zeta$ using purely geometrical considerations. Using the analytical expressions of the local strain rate $\dot{\epsilon}(x)$, we obtain expressions for the maximum (3D) or characteristic (2D) strain rate in a given unit cell and equate this with our definition of the nominal strain rate [Eq. (S11)]. Solving for $\zeta$ produces an expression that can be evaluated using known numerical constants and the medium porosity. Because the expressions for the local strain rate in the 2D arrays are unbounded [Figs. S3(e, f)], we instead define the characteristic strain rate as $\dot{\epsilon} = \frac{V_{\max}-V_{\min}}{L_{\min}-L_{\max}}$. For the 3D arrays, we determine expressions for the maximum local strain rate (marked in Figs. S4(e, f)) which produces an upper bound estimate for $\zeta$. Table S1 summarizes the theoretical $\zeta$ predictions. We find reasonable agreement with the experimentally-determined values for all four geometries. In summary, while $\zeta$ is a fitting parameter in our extensional resistance-modified power balance model, we show that the values of $\zeta$ for the ordered geometries have a physical basis by rationalizing the effective length scales for extensional strain in each unit cell geometry and comparing with theoretical predictions.

Derivation of modified Darcy's law incorporating an extensional resistance-per-stagnation point in ordered geometries

The analysis in § suggests that the extensional flow in the 2D pillar arrays and 3D sphere packings is predominantly controlled by stagnation points rather than contractions in the pore space. Accordingly, we seek to develop an alternative resistance model that captures the enhanced extensional flow resistance on a resistance-per-stagnation point basis.

We treat the optimized cross slot extensional rheometer (OSCER) [41, 42] as a model pore containing one free stagnation point (yellow $x$ marked in fig:S5a) at the center of the device. In planar extensional flow mode, flow is driven through two opposing inlets and exits through two perpendicular opposing outlets, thereby generating a strong extensional velocity gradient in the outflow channels. We measure the total pressure drop across the device, $\Delta P_{\mathrm{tot}}$, while applying a constant flow rate ramp. The device can also operate in a shear flow mode: one inlet and one outlet are closed off such that the fluid passes through a predominantly shearing corner flow. We measure the pressure drop associated with shearing in this mode as $\Delta P_{\mathrm{sh}}$. Assuming the excess pressure drop across the device results from purely extensional effects, we then obtain the extensional pressure drop in the device as $\Delta P_{\mathrm{ex}} = \Delta P_{\mathrm{tot}} - \Delta P_{\mathrm{sh}}$. Figure S5(b) plots the ratio of the extensional to shear pressure drop in the OSCER device as a function of the imposed extensional Weissenberg number, $\mathrm{Wi}_{E}$, defined as the product of the longest polymer relaxation time ($\lambda$) and a nominal strain rate in the device ($\dot{\epsilon}_{\mathrm{nom}}$). Error bars represent one standard deviation from time-averaging the pressure drop signal. The shaded gray region indicates the onset of an elastic instability characterized by asymmetry in the inflow/outflow velocities and a mobile stagnation point [43, 44, 45, 46]. As $\mathrm{Wi}_{E}$ increases, the extensional contribution to the pressure drop increases greatly, resulting from the onset of polymer stretching by the extensional velocity gradient originating at the stagnation point. We fit a power law to the pressure ratio measurements, $\frac{\Delta P_{\mathrm{ex}}}{\Delta P_{\mathrm{sh}}} \times100\% = 10.4\mathrm{Wi}_{E}^{1.8}$, which we take to model the extensional flow resistance-per-stagnation point. Note that the specific power law fitting parameters are fluid-dependent.

We now return to Darcy’s law and explicitly separate the pressure drop into extensional and shear components,

Factoring out the shear pressure drop, rearranging for the apparent viscosity, and normalizing by the bulk shear viscosity produces,

The first bracketed term is equivalent to the power balance model accounting for the unstable viscous dissipation only (4), and the right bracketed term is a factor accounting for the enhanced extensional resistance from one stagnation point. To correctly account for multiple stagnation points, we multiply the pressure ratio by the number of stagnation points per unit cell, $n_{\mathrm{SP}}$, which is known a priori for the ordered geometries. The resulting modified Darcy’s law accounting for extensional flow resistance is,

This model can be evaluated directly using the local flow fields ($\langle\chi\rangle_{t,V}$), fluid rheology ($\eta_{\mathrm{I}}$), and measurement of the pressure ratio in the OSCER device $\left(\frac{\Delta P_{\mathrm{ex}}}{\Delta P_{\mathrm{sh}}}\right)$ without requiring any additional fitting parameters given that $n_{\mathrm{SP}}$ is known. In the 2D hexagonal arrays, $n_{\mathrm{SP}}=6$. Table S2 provides a breakdown of the number of stagnation points per unit cell in the 3D packings, accounting for upstream and downstream stagnation points located at both the consolidations and the grain surfaces.

Geometry	Consolidations	Grains	Total
SC (3D)	8	0	8
BC (3D)	12	10	22

Disordered packings: power balances

Disordered sphere packings

As many porous media in nature and in industrial applications are inherently disordered, we now investigate the applicability of the extensional resistance-modified power balance models towards predicting the flow resistance in two types of disordered porous media: lightly-sintered disordered bead packings Figure S6(a) and lightly-sintered crushed glass packings Figure S7(a).

Following our previous work [2, 47, 48], we fabricate disordered bead packings by lightly sintering randomly-packed borosilicate glass beads ($D_{\mathrm{p}}=300$–$355\,\mu\mathrm{m}$; Mo-Sci) in a $3\times3$ mm rectangular quartz capillary. Figure S6(a) presents a sample confocal micrograph of a $z$-slice within one packing showing the disordered nature of the pore space geometry. Table S3 summarizes the relevant medium parameters. We repeat the same experimental flow protocol described in §, with one difference: due to the disordered nature of the pore space—where there are no longer characteristic unit cells—we randomly select and image 15 different pores within the porous medium in order to more accurately capture the pore-scale flow behavior.

Medium	$\phi$	$k$ ($\mu\mathrm{m}^{2}$)	$L$ (cm)	$A$ ($\mathrm{mm}^{2}$)
1	0.321	$82.9 \pm0.8$	2.62	9
2	0.339	$93.6 \pm0.8$	2.4	9
3	0.332	$83.3 \pm0.7$	2.02	9

Figures S6(b-d) show the extensional resistance-modified power balances (7) for 3 different experimental replicates in the disordered sphere packings, presented in order of decreasing model accuracy. The top row of plots corresponds to the local apparent extensional viscosity power balance model using a locally-weighted Trouton ratio (7). The power balance captures the macroscopic flow resistance in medium 1; in media 2 and 3, the model–data discrepancy is larger. The kink in the model curves for media 2 and 3 reflects mis-prediction of the onset of the extensional-resistance term: the abrupt uptick at $\mathrm{Wi}\approx3.5$ (medium 2) and 4.5 (medium 3) signals the crossover from a regime dominated by unstable viscous dissipation to one dominated by the extensional-resistance term. The relative weights of the two resistance terms are set by the critical Weissenberg number ($\mathrm{Wi}_{\mathrm{c}}$) for the onset of elastic instability and by the flow-field-weighted Trouton ratio ($\sum\mathrm{Tr}\,p(\mathrm{Tr})\,\delta\mathrm{Tr}$). Because ($\mathrm{Wi}_{\mathrm{c}}$) can be determined reliably from flow-field fluctuations, we attribute the model variability primarily to uncertainty in the weighted apparent Trouton ratio.

Consistent with prior measurements by Browne and Datta [2], we find that the unstable viscous dissipation power balance (black curves) captures the macroscopic flow resistance near the onset of flow thickening but progressively underpredicts flow thickening at higher $\mathrm{Wi}$—motivating the inclusion of an extensional-resistance term, as in the ordered geometries. Notably, the extent of underprediction varies across the three media (e.g., the unstable-dissipation balance tracks the data more closely near onset for medium 3, more similar to the results of Browne and Datta [2], than for medium 1, which differs more substantially). Two additional features of our data point to medium-to-medium heterogeneity as the origin of this variability: (i) the large-$\mathrm{Wi}$ plateau values of $\eta_{\mathrm{app}}/\eta_{\mathrm{I}}$ differ across replicates by more than experimental noise alone, and (ii) the $\mathrm{Wi}$ at which flow-thickening onsets shifts modestly between media prepared with the same polymer solution. We attribute this variability to the “porous individualism” of disordered packings [2]: each finite-volume realization samples a different local distribution of stagnation points, Hencky strains, and chaotic-flow regions. By contrast, ordered geometries yield reproducible flow behavior across unit-cell locations sufficiently far from the inlet and outlet, because of the periodic unit-cell geometry. We therefore expect that a statistically robust evaluation of the power-balance models in disordered media will require larger pore-scale samples ($n \gg15$ pores) than used here.

The fitted $\zeta$ values are similar across the 3 disordered packings ($\zeta= 1.28 - 1.57$), corresponding to a characteristic extensional length scale $\sim0.6 - 0.8\ D_p$. This length scale falls between the pore throat diameter ($0.16D_p$ [49, 50]) and grain diameter ($D_p$), which correspond to extensional flow generated in pore throat constrictions and stagnation points arising at grain surfaces, respectively.

The bottom row in Figure S6 (b-d) presents the resistance-per-stagnation point power balance model (14) for the 3 disordered packings. Here, we simply use $n_{\mathrm{SP}}$ as a fitting parameter to estimate extensional effects as it is challenging to (1) define a representative unit cell volume in the disordered packings and (2) accurately identify the locations of stagnation points within the 3D pore space. The models appear to suitably capture the correct order of magnitude of the flow resistance; however, the development of this resistance model was specifically founded on the premise of primarily stagnation point-dominated extensional flow in the ordered geometries, which is in general not applicable to disordered packings with highly varying flow topology [36, 51, 52]. The fitted values of $n_{\mathrm{SP}}$ (number of stagnation points per characteristic unit cell) also vary greatly across different experiments and have seemingly large values compared to the ordered porous medium geometries.

Medium	$\phi$	$k$ ($\mu\mathrm{m}^{2}$)	$L$ (cm)	$A$ ($\mathrm{mm}^{2}$)
1	0.297	$365 \pm4$	2.91	9
2	0.321	$218 \pm2$	3.11	9
3	0.292	$233 \pm2$	2.94	9

Crushed glass packings: power balances

We next examine the polymer solution flow resistance in packings of crushed glass beads, which produce pore space morphologies more representative of geological porous media encountered in subsurface flow applications. Briefly, we manually crush 1 mm borosilicate beads, sieve the particulates to achieve a target size range of $300$-$400\,\mu\mathrm{m}$, and lightly sinter packings of the glass fragments in a $3 \times3$ rectangular quartz capillary. Figure S7(a) shows a sample confocal micrograph of a $z$-slice within one of the crushed glass packings—note the enhanced pore space heterogeneity compared to the disordered bead packings [Fig. S6(a)]. Table S4 summarizes the relevant medium parameters for three different experimental replicates. Figures S7(b-d) present the resulting extensional-modified power balance models for the crushed glass packings, where the top row shows the local apparent extensional viscosity model [Eq. (S5)] and the bottom row shows the resistance-per-stagnation point model [Eq. (S12)], using $n_{\mathrm{SP}}$ as a fitting parameter. Experiments in panels (b-d) for different media are presented in order of decreasing accuracy. For media 1 and 2, the local apparent extensional viscosity model recapitulates the flow resistance in the crushed packings. The kink in the model curve for medium 2 results from the weighted Trouton ratio calculation, where the data for interpolation are not entirely smooth due to experimental variability. In contrast to media 1 and 2, the flow resistance of medium 3 is not captured well by the local extensional viscosity model, as there seems to be a mismatch in the onset of the extensional resistance term predicted by the model (see discussion in §). The values of the fitting parameter $\zeta$ in the crushed packings are more varied compared to the disordered bead packings, which may reflect the enhanced heterogeneity of the pore space and subsequently a broadening in the distribution of length scales for extensional flow.

As discussed in §, applying the resistance-per-stagnation point model to disordered packing geometries is no longer appropriate as the flow topology in the crushed packings likely consists of a complex mix of shear and extensional flow—even from pore-to-pore at the same imposed flow rate. We present the resistance-per-stagnation point models here for thoroughness in the bottom row of Figure S7(b-d), using $n_{\mathrm{SP}}$ as a fitting parameter rather than as a known geometrical parameter. The resulting models again show large variability in accuracy, and the fitted values of $n_{\mathrm{SP}}$ also vary greatly across different experiments.

In conclusion, our analyses in the disordered bead packings and crushed glass packings highlight the increased complexity in translating pore-scale flow behavior to the macroscopic flow resistance due to (1) increased spatial heterogeneity at the pore-scale requiring adequate flow sampling and (2) the mixed flow topology, where the flow is not dominated by extensional flow as we have found for the ordered porous medium geometries. Nevertheless, we expect that multiple mechanisms of flow resistance still contribute to the flow-thickening phenomenon in disordered porous media and expanding on the proposed resistance models to improve our prediction of these contributions will be a subject of future work.

Counting stagnation points in disordered packings

The fitted $n_{\mathrm{SP}}$ values for the disordered packings are an order of magnitude or more larger than the counted $n_{\mathrm{SP}}$ for the ordered geometries: $n_{\mathrm{SP}} = 280, 182, 108$ across the three disordered bead packings and $n_{\mathrm{SP}} = 37, 12, 41$ across the three crushed-glass packings. Three observations help interpret these values. First, in the absence of a well-defined unit cell, $n_{\mathrm{SP}}$ no longer represents a counted geometric integer but an effective parameter that combines the stagnation-point density per unit volume with the per-point contribution to the macroscopic resistance. Second, disordered packings contain many more stagnation points per unit volume than the ordered geometries studied here—grain-grain contacts occur throughout the medium, and the tortuous, irregular flow paths generate stagnation regions wherever streamlines bend sharply—making the order-of-magnitude increase in fitted $n_{\mathrm{SP}}$ relative to the ordered case geometrically plausible. Third, the fitted $n_{\mathrm{SP}}$ is smaller for the crushed-glass packings than for the bead packings, despite the sharper edges and more irregular geometry of the former; the contrast likely reflects differences in pore-scale geometry, permeability, and porosity between the two media, though a quantitative decomposition of these effects into a stagnation-point density and a per-point efficacy is beyond the scope of this work. The replicate-to-replicate variability within each medium type (a factor of approximately 3 in the bead packings, larger in the crushed-glass packings) reflects the “porous individualism” of disordered packings: each finite-volume realization samples a different local distribution of stagnation points and Hencky strains. (2) from the main text thus applies in disordered media only with $n_{\mathrm{SP}}$ as a fit parameter and with appropriate caution, in contrast to its parameter-free transferability across ordered geometries. Flows in disordered porous media are a complex mixture of shear and extension [52], which violates the assumptions underlying (2) of the main text.

Resistance from shear-generated normal stresses

Distinguishing between shear and extensionally-generated polymeric normal stresses

In (6), the term $(\tau_{xx}-\tau_{yy})\frac{\partial u}{\partial x}$ contains the tensile stress difference $\tau_{xx}-\tau_{yy}$. In our derivation, we substitute this stress difference with the apparent extensional viscosity arising from extensional flow. However, these tensile stresses may also arise from shearing flow, described by the first normal stress difference, $N_1$. This distinction was delineated by James et al. [39, 53]. Accordingly, we can introduce resistance contributions from both $\eta_{\mathrm{E}} = (\tau_{xx}-\tau_{yy})/\dot{\epsilon}$ and $N_1 = (\tau_{xx}-\tau_{yy})_{\mathrm{shear}}$ by separating out the stresses generated by shear and extensional flow. The revised resistance model follows as,

A first approximation to evaluating the normal stress difference term in (15) is simply to evaluate the steady-state $N_1$ power law model obtained in bulk rheology using the instantaneous local pore-scale shear rate fields. Additionally, the values of $\frac{\partial u}{\partial x}$ can take on both negative and positive signs: negative indicates compression and positive indicates extension. For positive $N_1$, the sign of the product $N_1\frac{\partial u}{\partial x}$ then determines whether the shear-generated tensile stress dissipates or absorbs work from the flow. Figure S8(a-d) show the average $E = N_1\frac{\partial u}{\partial x}$ values as a function of Wi measured from the flow fields in the different ordered geometries, where negative values indicate net energy dissipation, while positive values indicate net energy accumulation due to storage of elastic stresses in the polymer chains. For all medium geometries, the net dissipation increases with Wi; however, the magnitude of this dissipation is quite small.

Here, we estimate the maximum possible upper bound to the dissipation from $N_1$ by taking the absolute value of the measured gradient $\frac{\partial u}{\partial x}$. Figure S8(e-h) present the resulting $E = N_1\left|\frac{\partial u}{\partial x}\right|$ curves as a function of Wi. We use these upper bound estimates to evaluate the full resistance model in (15). Figure S8(i-l) show the ratio of the resistances contributed by the normal stress difference (fourth term on the right in (15)) and extensional viscosity (third term on the right in Eq. (S13)). At low Wi, the $N_1$ contribution dominates. As Wi increases, this ratio rapidly becomes dominated by the extensional viscosity resistance across all of the ordered geometries. Note that the data shown correspond to an upper bound estimate of $N_1$; hence, the true contribution from $N_1$ is likely lower. We thus conclude that, in the ordered porous medium geometries, the tensile stresses arising from extensional flow govern the flow-thickening behavior compared to shear-generated tensile stresses.

Assessment of transient versus steady estimates of $N_1$

The first normal stress difference is fundamentally a transient property that depends on the history of shearing. Thus, evaluating the steady-state $N_1$ power law model from bulk rheology using the local instantaneous pore-scale shear rate fields may be inaccurate. To capture the transient nature of $N_1$ in evaluation of the resistance model, we first perform transient bulk shear rheology measurements by imposing a step shear strain and recording the stress response over time at various shear rates (Figure S9(a)). The range of resolvable shear rates is limited on the lower end by the instrument resolution and on the higher end by the onset of an elastic instability in the cone-plate geometry. Rescaling the $N_1$ data by the steady-state plateau value $N_{1,p}$ and plotting against the accumulated shear strain $\gamma= \dot{\gamma}t$, we observe a collapse in the transient $N_1$ response in Figure S9(b). This collapse can be used as a mapping from measured shear strain to transient $N_1$, where the shear strains can be measured from the spatiotemporal flow fields as $\gamma= \int_{t=0}^{t=\lambda} \dot{\gamma}\,\mathrm{d}t$.

Figure S9(c-f) present distributions of the measured shear strains, $\gamma$, and extensional Hencky strains, $\epsilon_{\mathrm{Hencky}}$, in each of the four ordered porous medium geometries. Generally, both types of strains increase with Wi. Importantly, comparing the maximum values of $\gamma$ with the transient rheology in Figure S9(b), we find that the accumulated shear strains in the pore-scale flow fields ($\gamma\approx10 - 20$) are smaller than the nominal shear strain ($\gamma\approx300$) for $N_1$ to reach a steady state. This result suggests that the calculations in Figure S8 are indeed an upper bound estimate for the resistance contributed by shear-generated normal stresses, neglecting the stress overshoot behavior at around 100 strain units. Further, while the accumulated strain distributions indicate that accumulated shear strains are larger in magnitude than the extensional strains, the full power balance calculation presented in Figure S8(i-l) indicate that energy dissipation by the associated transient extensional stresses is larger than the shear stresses, particularly with increasing Wi. We thus conclude that—in the ordered porous medium geometries—the polymeric tensile stresses $\tau_{xx} - \tau_{yy}$ primarily arise from extensional flow rather than shearing flow. We expect that this may not be the case in other types of porous medium geometries, as altogether our results indicate that the flow topology is intricately linked to the pore space geometry and thus the macroscopic flow resistance.

References