1. Introduction
Many key environmental, industrial and energy processes – such as remediation of contaminated groundwater aquifers (Smith et al. Reference Smith, Silva, Munakata-Marr and McCray2008; Hartmann et al. Reference Hartmann2021), recovery of oil from subsurface reservoirs (Durst, Haas & Kaczmar Reference Durst, Haas and Kaczmar1981; Sorbie Reference Sorbie2013) and extraction of heat from geothermal reservoirs (Di Dato et al. Reference Di Dato, D'Angelo, Casasso and Zarlenga2021) – rely on the injection of a fluid into a subsurface porous medium. Such media are formed by sedimentary processes, often leading to vertically layered strata of distinct pore sizes oriented along the direction of macroscopic flow (Freeze Reference Freeze1975; Dagan Reference Dagan2012). The permeability differences between these strata cause uneven fluid partitioning across them, with preferential flow through higher-permeability regions and ‘bypassing’ of lower-permeability regions (Lake & Hirasaki Reference Lake and Hirasaki1981; Di Dato et al. Reference Di Dato, D'Angelo, Casasso and Zarlenga2021). This flow heterogeneity reduces the efficacy of contaminant remediation, oil recovery and heat extraction from bypassed regions – necessitating the development of new ways to spatially homogenize the flow.
Polymer additives have a long history of use in such applications to increase the injected fluid viscosity and thereby suppress instabilities, like viscous fingering, at immiscible (e.g. water–oil) interfaces (Durst et al. Reference Durst, Haas and Kaczmar1981; Smith et al. Reference Smith, Silva, Munakata-Marr and McCray2008; Sorbie Reference Sorbie2013). However, for typically used polymer solutions that do not have appreciable elasticity, this process of conformance control still suffers from the issue of uneven partitioning of flow across different strata due to differences in permeability. Quantitatively, the superficial velocity in a given stratum 
$i$ is given by Darcy's law, representing each stratum as a homogeneous medium with uniformly disordered pores of a single mean size: 
$U_{i}\equiv Q_{i}/A_{i}=(\Delta P/L)k_{i}/\eta _{{app}}$, where 
$Q_{i}$ is the volumetric flow rate through the stratum, 
$\Delta P$ is the pressure drop across a length 
$L$ of the parallel strata, 
$A_{i}$ and 
$k_{i}$ are the cross-sectional area and permeability of the stratum, respectively, and 
$\eta _{{app}}$ is the ‘apparent viscosity’ of the polymer solution quantifying the macroscopic resistance to flow through the tortuous pore space. For non-elastic polymer solutions, 
$\eta _{{app}}$ is simply given by the dynamic shear viscosity 
$\eta$ of the solution, and is typically not strongly dependent on flow rate or porous medium geometry. Therefore, differences in 
$k_{i}$ result in differences in 
$U_{i}$ between strata – leading to uneven partitioning of the flow across the entire stratified medium.
Conversely, the apparent viscosity of highly elastic polymer solutions (e.g. with molecular weights 
${\gtrsim }1\,\text {MDa}$) can depend on flow rate. For many such solutions, 
$\eta _{{app}}$ strongly increases above a threshold flow rate in a homogeneous porous medium, even though 
$\eta$ of the bulk solution decreases with increasing shear rate (Marshall & Metzner Reference Marshall and Metzner1967; James & McLaren Reference James and McLaren1975; Chauveteau & Moan Reference Chauveteau and Moan1981; Durst & Haas Reference Durst and Haas1981; Durst et al. Reference Durst, Haas and Kaczmar1981; Kauser et al. Reference Kauser, Dos Santos, Delgado, Muller and Saez1999; Haward & Odell Reference Haward and Odell2003; Odell & Haward Reference Odell and Haward2006; Zamani et al. Reference Zamani, Bondino, Kaufmann and Skauge2015; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016; Skauge et al. Reference Skauge, Zamani, Gausdal Jacobsen, Shaker Shiran, Al-Shakry and Skauge2018; Ibezim, Poole & Dennis Reference Ibezim, Poole and Dennis2021). Direct visualization of the flow in a homogeneous medium (Browne & Datta Reference Browne and Datta2021) recently established that this anomalous increase coincides with the onset of a purely elastic flow instability arising from the buildup of polymer elastic stresses during transport (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990; McKinley, Pakdel & Öztekin Reference McKinley, Pakdel and Öztekin1996; Pakdel & McKinley Reference Pakdel and McKinley1996; Shaqfeh Reference Shaqfeh1996; Burghelea et al. Reference Burghelea, Segre, Bar-Joseph, Groisman and Steinberg2004; Rodd et al. Reference Rodd, Cooper-White, Boger and McKinley2007; Afonso, Alves & Pinho Reference Afonso, Alves and Pinho2010; Galindo-Rosales et al. Reference Galindo-Rosales, Campo-Deaño, Pinho, Van Bokhorst, Hamersma, Oliveira and Alves2012; Zilz et al. Reference Zilz, Poole, Alves, Bartolo, Levaché and Lindner2012; Ribeiro et al. Reference Ribeiro, Coelho, Pinho and Alves2014; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016; Machado et al. Reference Machado, Bodiguel, Beaumont, Clisson and Colin2016; Kawale et al. Reference Kawale, Marques, Zitha, Kreutzer, Rossen and Boukany2017; Sousa, Pinho & Alves Reference Sousa, Pinho and Alves2018; Browne, Shih & Datta Reference Browne, Shih and Datta2019; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019a; Browne, Shih & Datta Reference Browne, Shih and Datta2020; Walkama, Waisbord & Guasto Reference Walkama, Waisbord and Guasto2020; Haward, Hopkins & Shen Reference Haward, Hopkins and Shen2021). Such an instability gives rise to random flow fluctuations that are, in some cases, reminiscent of inertial turbulence despite the vanishingly small Reynolds numbers 
${Re}$ (Groisman & Steinberg Reference Groisman and Steinberg2000; Pan et al. Reference Pan, Morozov, Wagner and Arratia2013; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019a; Datta et al. Reference Datta2022) – contributing added viscous dissipation that, at least in some cases, primarily generates this anomalous increase in 
$\eta _{{app}}$ (Browne & Datta Reference Browne and Datta2021). In a stratified medium, this flow rate dependence of 
$\eta _{{app},i}$ in each stratum may provide an avenue to break the proportionality between 
$k_{i}$ and 
$U_{i}$, potentially mitigating the uneven partitioning of the flow across strata. However, this possibility remains unexplored; indeed, it is still unknown how exactly this elastic instability arises in each stratum.
Here, we demonstrate that this elastic flow instability can help homogenize flow in stratified porous media. Using pore-scale confocal microscopy and macro-scale imaging of passive scalar transport, we visualize the flow in a model porous medium with two distinct parallel strata, imposing a constant flow rate 
$Q$ through the entire medium. For small 
$Q$, the flow in both strata is laminar, leading to the typical uneven partitioning of flow across the strata. Strikingly, for 
$Q$ above a threshold value, the instability arises solely in the higher-permeability stratum and fluid is redirected to the lower-permeability stratum, helping to homogenize the flow. Above an even larger threshold flow rate, the instability also arises in this lower-permeability stratum, suppressing this flow redirection – leading to a window of flow rates at which this homogenization arises. Guided by these findings, we develop a parallel-resistor model that treats each stratum 
$i$ as a homogeneous medium with specified 
$A_i$, 
$k_i$ and, therefore, 
$\eta _{{app},i}$, all coupled at the inlet and outlet. This model quantitatively captures the overall pressure drop across the stratified medium as well as the observed flow redirection with varying flow rate. It also elucidates the underlying cause of this redirection. In particular, above the first threshold flow rate, preferential flow causes the elastic instability to arise solely in the higher-permeability stratum. The corresponding increase in the resistance to flow, as quantified by 
$\eta _{{app},i}$, redirects flow towards the lower-permeability stratum. Above the larger second threshold flow rate, the onset of the instability and corresponding increase in 
$\eta _{{app},i}$ in the lower-permeability stratum redirects flow back towards the higher-permeability stratum – yielding the experimentally observed optimum in flow homogenization. Finally, we generalize this model, establishing the operating conditions at which this homogenization is optimized for porous media with arbitrarily many strata. Thus, our work provides a new approach to homogenizing fluid and passive scalar transport in heterogeneous porous media. Since many naturally occurring media are stratified, we anticipate these findings to be broadly useful in environmental, industrial and energy processes.
2. Materials and methods
To investigate the spatial distribution of flow in a stratified porous medium, we use imaging at two different length scales: macro-scale (
${\sim }100$ s of pores, figure 1a) and pore-scale (
$\sim$1 pore, figure 1b). For clarity regarding nomenclature, we note that different forms of elastic flow instability may arise in different geometries (e.g. those with curvilinear vs parallel streamlines), possibly resulting in flows with distinct power spectra of flow velocity fluctuations (Groisman & Steinberg Reference Groisman and Steinberg2000; Fouxon & Lebedev Reference Fouxon and Lebedev2003; Steinberg Reference Steinberg2021, Reference Steinberg2022; Datta et al. Reference Datta2022). Clarifying the different forms of elastic instability that can arise, and the physics underlying the transition to and spectral features of each instability, is an interesting question for future research. However, the work described in this manuscript instead focuses on the flow resistance, and corresponding partitioning of flow across strata, arising from an elastic instability in porous media – not on the exact nature of the instability or the transition to unstable flow. We therefore use the term ‘elastic flow instability’ to refer more generally to polymer elasticity-generated flow instabilities at low 
$Re \ll 1$, independent of the specific details of the exact nature of the instability, but instead dependent on the flow thickening it causes as shown by our experimental measurements.
Figure 1. Macro-scale and pore-scale characterization of the flow of a highly elastic polymer solution in stratified porous media. (a) Schematic of our Hele-Shaw assembly, with two parallel strata made of close-packed glass beads of distinct sizes. We characterize the macro-scale flow by visualizing the transport of an injected dye, which acts as a passive scalar. (b) Schematic of our microfluidic assembly, again with two parallel strata made of close-packed glass beads of the same distinct sizes as in the Hele-Shaw assembly. We characterize the flow through direct pore-scale visualization of fluorescent tracer particle transport using confocal microscopy, combined with measurements of the pressure drop across the entire medium. (c) Shear stress 
$\sigma$ varies nearly linearly with the shear rate 
$\dot {\gamma }$, indicating that the solution approximates a Boger fluid; the dashed line shows the power-law fit 
$\sigma =A_s\dot {\gamma }^{\alpha _s}$ with 
$A_s=0.3428\pm 0.0002\,\mathrm {Pa}\,\mathrm {s}^{\alpha _s}$ and 
$\alpha _s=0.931\pm 0.001$. (d) First normal stress difference 
$N_1$ also increases with increasing shear rate 
$\dot {\gamma }$; the dashed line shows the power-law fit 
$N_1=A_n\dot {\gamma }^{\alpha _n}$ with 
$A_n=1.16\pm 0.03\,\mathrm {Pa}\,\mathrm {s}^{\alpha _n}$ and 
$\alpha _n=1.25\pm 0.02$. Error bars represent one standard deviation of three measurements.
2.1. Macro-scale experiments in a Hele-Shaw assembly
To characterize the macro-scale partitioning of flow, we fabricate an unconsolidated stratified porous medium in a Hele-Shaw assembly (figure 1a). We 3-D print an open-faced rectangular cell with spanwise (
$y$–
$z$-direction) cross-sectional area 
$A=3\,\text {cm}\times 0.4\,\text {cm}$ and streamwise (
$x$-direction) length 
$L=5\,\text {cm}$ using a clear methacrylate-based resin (FLGPCL04, Formlabs Form3). To ensure an even distribution of flow at the boundaries, both the inlet and outlet are split into three equally spaced streams 
${\approx }2\,\text {cm}$ from the entrance and exit of the porous medium, respectively, with the centre stream spanning the streamwise interface between the parallel strata. We then fill the cell with spherical borosilicate glass beads of distinct diameters arranged in parallel strata using a temporary partition, with bead diameters 
$d_p=1000$ to 1400 
$\mathrm {\mu }$m (Sigma Aldrich) and 212 to 255 
$\mathrm {\mu }$m (Mo-Sci) for the higher-permeability coarse (subscript 
$C$) and lower-permeability fine (subscript 
$F$) strata, respectively. The strata have equal cross-sectional areas 
$A_C\approx A_F\approx A/2$ and thus their area ratio 
$\tilde {A}\equiv A_C/A_F\approx 1$. Steel mesh with a 150 
$\mathrm {\mu }$m pore size cutoff placed over the inlet and outlet tubing prevents the beads from exiting the cell. We tamp down the beads for 30 min to form a dense random packing with a porosity 
$\phi _V\sim 0.4$ (Onoda & Liniger Reference Onoda and Liniger1990). We then screw the whole assembly shut with an overlying acrylic sheet cut to size, sandwiching a thin sheet of polydimethylsiloxane to provide a watertight seal. All experiments using this assembly are conducted at 
${\approx }21\,^{\circ }\mathrm {C}$.
For all macro-scale experiments, we use a Harvard Apparatus PHD 2000 syringe pump to first introduce the test fluid – either the polymer solution or the polymer-free solvent, which acts as a Newtonian control – at a constant flow rate 
$Q$ for at least the duration needed to fill the entire pore space volume 
$t_{PV}\equiv \phi _V AL/Q$ before imaging to ensure an equilibrated starting condition. We then visualize the macro-scale scalar transport by the fluid by introducing a step change in the concentration of a dilute dye (0.1 wt.% green food colouring, McCormick) and record the infiltration of the dye front using a DSLR camera (Sony 
$\alpha$6300), as shown in figure 1(a). To track the progression of the dye as it is advected by the flow, we determine the ‘breakthrough’ curve half-way along the length of the medium (
$x=L/2$) by measuring the dye intensity 
$C$ averaged across the entire medium cross-section, normalized by the difference in intensities of the final dye-saturated and initial dye-free medium, 
$C_f$ and 
$C_0$, respectively: 
$\tilde {C}\equiv (\langle C\rangle _y-\langle C_0\rangle _y)/(\langle C_f\rangle _y-\langle C_0\rangle _y)$ (figure 2b). This procedure averages out slight variations in the dye front along the 
$y$- and 
$z$-directions, which inevitably arise due to grain-scale spatial fluctuations in the pore geometry (Datta et al. Reference Datta, Chiang, Ramakrishnan and Weitz2013). For all breakthrough curves thereby measured, time 
$t$ is normalized using the time taken to reach this half-way point, 
$\tilde {t}\equiv t/(0.5t_{PV})$. Repeating this procedure for individual strata (subscript 
$i$) and tracking the variation of the streamwise position 
$X_i$ at which 
$\tilde {C}_i=0.5$ with time provides a measure of the superficial velocity 
$U_i=\mathrm {d}X_i/\mathrm {d}t$ in each stratum. In between tests at different flow rates, we flush the assembly with the dye-free solution for at least ten pore volumes to remove any residual dye.
Figure 2. Imaging reveals that an elastic polymer solution homogenizes the uneven flow between strata, coincident with the onset of the elastic flow instability in the coarser stratum. (a) Visualization of passive scalar transport by the polymer solution in a stratified Hele-Shaw assembly using a green dye. All images are taken at the same 
$\tilde {t}\equiv t/(0.5t_{PV})=2.5$, where time 
$t$ has been normalized by the time to fill half of the entire pore space volume. Due to the higher permeability of the coarse stratum (bottom), dye infiltrates faster than in the fine stratum (top). However, at the intermediate 
${Wi}_I=2.7$, this uneven partitioning of the flow is reduced. (b) Scalar breakthrough curves obtained by measuring the normalized dye concentration 
$\tilde {C}$ at the midpoint 
$x=L/2$ over time. Uneven flow partitioning at 
${Wi}_I=1.4$ leads to distinct jumps and prolongs 
$\tilde {C}$ to long times; by contrast, redirection of flow to the fine stratum at the intermediate 
${Wi}_I=2.7$ leads to more uniform and rapid breakthrough, shown by the smoother and earlier rise in 
$\tilde {C}(\tilde {t})$. This homogenization is mitigated at the even larger 
${Wi}_I=3.3$. (c) Streamline images of representative pores in a stratified microfluidic assembly; black circles are sections through the beads making up the solid matrix, white lines are time projections of the tracer particle pathlines that closely approximate the instantaneous flow streamlines. Imposed flow direction is from left to right. The flow homogenization at the intermediate 
${Wi}_I=2.7$ (c i,c iii) coincides with the onset of the elastic flow instability solely in the coarse stratum (bottom) – indicated by the emergence of spatio-temporal fluctuations in the flow, shown by the red overlay whose intensity is given by the standard deviation in pixel intensity over the course of the time series of images. The mitigation of this homogenization at the even larger 
${Wi}_I=3.3$ (c ii,c iv) coincides with the additional onset of the instability in the fine stratum, as well (top). (d) Map of the root mean square velocity fluctuations computed from PIV, which confirm that the flow becomes unstable first in the coarse stratum at the intermediate 
${Wi}_I=2.7$ (d i,d iii), and then also in the fine stratum at the even larger 
${Wi}_I=3.3$ (d ii,c iv). (e) Fraction of 10 randomly chosen pores observed in each stratum that exhibit unstable flow, defined as such by identifying whether fluid streamlines cross over the imaging duration. Only a small fraction of pores in the coarse stratum exhibit unstable flow at the intermediate 
${Wi}_I=2.7$, whereas a greater fraction of pores in both strata exhibit unstable flow at the larger 
${Wi}_I=3.3$ – corroborating the results shown in (c–d).
2.2. Pore-scale experiments in microfluidic assemblies
To gain insight into the pore-scale physics, we use experiments in consolidated microfluidic assemblies (figure 2b). We pack spherical borosilicate glass beads (Mo-Sci) in square quartz capillaries (
$A= 3.2\,{\rm mm} \times 3.2\,{\rm mm}$; Vitrocom), densify them by tapping and lightly sinter the beads – resulting in dense random packings again with 
$\phi _V\sim 0.4$ (Krummel et al. Reference Krummel, Datta, Münster and Weitz2013). We use this protocol to fabricate three different microfluidic media: a homogeneous higher-permeability coarse medium (
$d_p= 300$ to 
$355\,\mathrm {\mu }$m), a homogeneous lower-permeability fine medium (
$d_p= 125$ to 
$155\,\mathrm {\mu }$m) and a stratified medium with parallel higher-permeability coarse and lower-permeability fine strata, each composed of the same beads used to make the homogeneous media, again with equal cross-section areas, 
$\tilde {A}\approx 1$ (Datta & Weitz Reference Datta and Weitz2013; Lu et al. Reference Lu, Pahlavan, Browne, Amchin, Stone and Datta2020). We measure the fully developed pressure drop 
$\Delta P$ across each medium using an Omega PX26 differential pressure transducer. All experiments on the microfluidic assemblies are conducted at 
${\approx }21\,^{\circ }\mathrm {C}$.
For all pore-scale experiments, before each experiment, we infiltrate the medium to be studied first with isopropyl alcohol (IPA) to prevent trapping of air bubbles and then displace the IPA by flushing with water. We then displace the water with the miscible polymer solution, seeded with 5 ppm of fluorescent carboxylated polystyrene tracer particles (Invitrogen), 
$D_t=200\,\mathrm {nm}$ in diameter. This solution is injected into the medium at a constant volumetric flow rate 
$Q$ using Harvard Apparatus syringe pumps – a PHD 2000 for 
$Q>1\,\text {ml}\,\text {h}^{-1}$ or a Pico Elite for 
$Q<1\,\text {ml}\,\text {h}^{-1}$ – for at least 3 h to reach an equilibrated state before flow characterization. After each subsequent change in 
$Q$, the flow is given 1 h to equilibrate before imaging. We monitor the flow in individual pores using a Nikon A1R+ laser scanning confocal fluorescence microscope with a 488 nm excitation laser and a 500–550 nm sensor detector; the tracer particles have excitation between 480 and 510 nm with an excitation peak at 505 nm, and emission between 505 and 540 nm with an emission peak at 515 nm. These particles are faithful tracers of the underlying flow field since the Péclet number 
${Pe}\equiv (Q/A)D_t/\mathcal {D}>10^5\gg 1$, where 
$\mathcal {D}=k_BT/(3{\rm \pi} \eta D_t)=6\times 10^{-3}\,\mathrm {\mu }\mathrm {m}^2\,\mathrm {s}^{-1}$ is the Stokes–Einstein particle diffusivity. We then visualize the flow using a 10
$\times$ objective lens with the confocal resonant scanner, obtaining successive 8 
$\mathrm {\mu }\mathrm {m}$-thick optical slices at a 
$z$ depth hundreds of 
$\mathrm {\mu }\mathrm {m}$ within the medium. Our imaging probes an 
$x$–
$y$ field of view 
$159\,\mathrm {\mu }\mathrm {m}\times 159\,\mathrm {\mu }\mathrm {m}$ at 60 frames per second for pores with 
$d_p=125$ to 
$155\,\mathrm {\mu }\mathrm {m}$ or 
$318\,\mathrm {\mu }\mathrm {m}\times 318\,\mathrm {\mu }\mathrm {m}$ at 30 frames per second for pores with 
$d_p=300$ to 
$355\,\mathrm {\mu }\mathrm {m}$.
To monitor the flow in the different pores over time, we use an ‘intermittent’ imaging protocol. Specifically, we record the flow in multiple pores chosen randomly throughout each medium (19 and 20 pores of the homogeneous coarse and fine media, respectively) for 2 s long intervals every 4 min over the course of 1 h. For the experiments in homogeneous fine and stratified media, we also complement this protocol with ‘continuous’ imaging in which we monitor the flow successively in 10 pores of the homogeneous fine medium for 5 min long intervals each. For ease of visualization, we intensity average the successive images thereby obtained over a time scale 
$\approx 2.5\,\mathrm {\mu }\mathrm {m}/(Q/A)$ (figure 2c), producing movies of the tracer particle pathlines that closely approximate the instantaneous flow streamlines. We also measure the instantaneous 
$x$–
$y$ velocities 
$\boldsymbol {u}$ using particle image velocimetry (PIV) (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014). Subtracting off the temporal mean, indicated by 
$\langle \rangle _t$, in each pixel yields the velocity fluctuation 
$\boldsymbol {u}'=\boldsymbol {u}-\langle \boldsymbol {u}\rangle _t$; we then define an unstable region of the flow as one in which the root mean square velocity fluctuation 
$\sqrt {\langle |\boldsymbol {u}'|^2\rangle _t}/\langle |\boldsymbol {u}|\rangle _t>0.3$.
2.3. Permeability measurements
For each medium, we determine the permeability via Darcy's law using experiments with pure water. For the microfluidic assemblies, we obtain 
$k_C=79\,\mathrm {\mu }\mathrm {m}^2$ and 
$k_F=8.6\,\mathrm {\mu }\mathrm {m}^2$ for the homogeneous coarse and fine media, respectively – comparable to our previously measured values on similar media (Krummel et al. Reference Krummel, Datta, Münster and Weitz2013) and to the prediction of the established Kozeny–Carman relation (Philipse & Pathmamanoharan Reference Philipse and Pathmamanoharan1993). The permeability ratio between the two strata is then 
$\tilde {k}\equiv k_C/k_F\approx 9$. The measured permeability for the entire stratified porous medium is 
$k=32\,\mathrm {\mu }\mathrm {m}^2$, in reasonable agreement with the prediction obtained by considering the strata as separated homogeneous media providing parallel resistance to flow, 
$k\approx \tilde {A} k_C+(1-\tilde {A})k_F\approx 44\,\mathrm {\mu }\mathrm {m}^2$.
The permeability of an isolated stratum in a stratified medium varies as 
$\sim d_p^2$, similar to a homogeneous porous medium. Hence, for the Hele-Shaw assembly, we estimate the permeability of each stratum by scaling 
$k_C$ and 
$k_F$ with the differences in bead size. We thereby estimate 
$k\approx 440\,\mathrm {\mu }\mathrm {m}^2$ (
$\tilde {k}\approx 26$) for the entire stratified medium, in reasonable agreement with the measured 
$k=270\,\mathrm {\mu }\text {m}^2$.
For both assemblies, we define a characteristic shear rate of the entire medium 
$\dot {\gamma }_I\equiv Q/(A\sqrt {\phi _V k})$ as the ratio between the characteristic pore flow speed 
$Q/(\phi _V A)$ and length scale 
$\sqrt {k/\phi _V}$ (Zami-Pierre et al. Reference Zami-Pierre, De Loubens, Quintard and Davit2016; Berg & van Wunnik Reference Berg and van Wunnik2017). Our experiments explore the range 
$\dot {\gamma }_I\approx 0.2$ to 
$26\,\mathrm {s}^{-1}$.
2.4. Polymer solution rheology
The polymer solution approximates a Boger fluid composed of dilute 300 ppm 18 MDa partially hydrolysed polyacrylamide dissolved in a viscous aqueous solvent composed of 6 wt.% ultrapure milliPore water, 82.6 wt.% glycerol (Sigma Aldrich), 10.4 wt.% dimethylsulfoxide (Sigma Aldrich) and 1 wt.% NaCl. This solution is formulated to precisely match its refractive index to that of the glass beads – thus rendering each medium transparent when saturated. From intrinsic viscosity measurements the overlap concentration is 
$c^{*}\approx 0.77/[\eta ]=600\pm 300\,\mathrm {ppm}$ (Browne & Datta Reference Browne and Datta2021) and the radius of gyration is 
$R_g\approx 220\,\mathrm {nm}$ (Rubinstein et al. Reference Rubinstein and Colby2003) and therefore our experiments use a dilute polymer solution at 
$\approx 0.5$ times the overlap concentration. The shear stress 
$\sigma (\dot {\gamma }_I)=A_s\dot {\gamma }^{\alpha _s}$ and first normal stress difference 
$N_1(\dot {\gamma }_I)=A_n\dot {\gamma }^{\alpha _n}$ are measured in an Anton Paar MCR301 rheometer, using a 
$1^{\circ }$ 5 cm diameter conical geometry set at a 50 
$\mathrm {\mu }$m gap, yielding the best-fit power laws 
$A_s=0.3428\pm 0.0002\,\mathrm {Pa}\,\mathrm {s}^{\alpha _s}$ with 
$\alpha _s=0.931\pm 0.001$ and 
$A_n=1.16\pm 0.03\,\mathrm {Pa}\,\mathrm {s}^{\alpha _n}$ with 
$\alpha _n=1.25\pm 0.02$ (figure 1c,d). All rheological measurements are temperature controlled with a Peltier plate at 
$20.0\pm 0.1\,^{\circ }\mathrm {C}$.
These measurements then enable us to calculate a characteristic interstitial Weissenberg number, which characterizes the role of polymer elasticity in the flow by comparing the magnitude of elastic and viscous stresses, 
${Wi}_I\equiv N_1(\dot {\gamma }_I)/(2\sigma (\dot {\gamma }_I))$, as commonly defined (Pan et al. Reference Pan, Morozov, Wagner and Arratia2013; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019a,Reference Qin, Salipante, Hudson and Arratiab; Browne et al. Reference Browne, Shih and Datta2020; Browne & Datta Reference Browne and Datta2021; Ibezim et al. Reference Ibezim, Poole and Dennis2021; Datta et al. Reference Datta2022). In our experiments this quantity exceeds unity, ranging from 1 to 5.5, suggesting that elastic flow instabilities likely arise in the flow (Larson et al. Reference Larson, Shaqfeh and Muller1990; Pakdel & McKinley Reference Pakdel and McKinley1996; Shaqfeh Reference Shaqfeh1996; Rodd et al. Reference Rodd, Cooper-White, Boger and McKinley2007; Afonso et al. Reference Afonso, Alves and Pinho2010; Galindo-Rosales et al. Reference Galindo-Rosales, Campo-Deaño, Pinho, Van Bokhorst, Hamersma, Oliveira and Alves2012; Zilz et al. Reference Zilz, Poole, Alves, Bartolo, Levaché and Lindner2012; Pan et al. Reference Pan, Morozov, Wagner and Arratia2013; Ribeiro et al. Reference Ribeiro, Coelho, Pinho and Alves2014; Sousa et al. Reference Sousa, Pinho and Alves2018; Browne et al. Reference Browne, Shih and Datta2019, Reference Browne, Shih and Datta2020; Browne & Datta Reference Browne and Datta2021) – as we directly verify using flow visualization, detailed further below.
Polymer solutions such as that used here are typically characterized by a broad spectrum of relaxation times, which can be challenging to fully characterize experimentally (Perkins, Smith & Chu Reference Perkins, Smith and Chu1997; Liu, Jun & Steinberg Reference Liu, Jun and Steinberg2007, Reference Liu, Jun and Steinberg2009). For simplicity, we therefore follow Shaqfeh (Reference Shaqfeh1996) and describe our solution using a single polymer relaxation time 
$\lambda \equiv {\lim }_{\dot {\gamma }\rightarrow 0}~{\varPsi _1}/({2\eta _p})$, where 
$\varPsi _1\equiv N_1/\dot {\gamma }^2$ is the first normal stress coefficient, 
$\eta _p\equiv \eta -\eta _s$ is the polymer contribution to the solution viscosity and 
$\eta _s=226.8\pm 0.3\,\mathrm {mPa}\,\mathrm {s}$ is the viscosity of the polymer-free solvent. Using the lowest shear rates where each value is accessible on our rheometer, we estimate 
${\lim }_{\dot {\gamma }\rightarrow 0}\varPsi _1\approx 192\pm 7\,\mathrm {mPa}\,\mathrm {s}^2$ and 
${\lim }_{\dot {\gamma }\rightarrow 0}\eta \approx 427\pm 4\,\mathrm {mPa}\,\mathrm {s}$, yielding 
$\lambda \approx 480\pm 30\,\text {ms}$. This value of a characteristic relaxation time is in good agreement with previously reported relaxation times for similar polymer and solvent compositions (Groisman & Steinberg Reference Groisman and Steinberg2000; Galindo-Rosales et al. Reference Galindo-Rosales, Campo-Deaño, Pinho, Van Bokhorst, Hamersma, Oliveira and Alves2012; Pan et al. Reference Pan, Morozov, Wagner and Arratia2013; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016; Walkama et al. Reference Walkama, Waisbord and Guasto2020), although we expect the true longest relaxation time of the solution to be larger than this value.
Finally, we also characterize the role of inertia with the Reynolds number 
${Re}=\rho Ud_p/\eta (\dot {\gamma }_I)$, which quantifies the ratio of inertial to viscous stresses for a fluid with density 
$\rho$. In our experiments this quantity ranges from 
${Re}=2\times 10^{-7}$ to 
$2\times 10^{-5}\ll 1$, indicating that inertial effects are negligible.
3. Results
3.1. Polymer solution homogenizes flow above a threshold Weissenberg number, coinciding with the onset of the elastic flow instability
We use our stratified Hele-Shaw assembly to characterize the uneven partitioning of flow between strata at the macro-scale. First, we impose a small flow rate 
$Q=3\,{\rm ml}\,{\rm h}^{-1}$ corresponding to 
${Wi}_I=1.4$ – below the onset of the elastic flow instability at 
${Wi}_I\approx 2.6$ for homogeneous media (Browne & Datta Reference Browne and Datta2021). As is the case with Newtonian fluids, we observe preferential flow through the coarse stratum, indicated by the infiltrating dye front in figure 2(a i) and in supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.337. The infiltration of dye at different rates through the strata produces two distinct steps in the breakthrough curve (blue line in figure 2b): the first jump from 
$\tilde {C}\approx 0$ to 
$0.4$ from 
$0<\tilde {t}\lesssim 3$ corresponds to fluid infiltration of the coarse stratum, and the second jump from 
$\tilde {C}\approx 0.4$ to 
$0.8$ from 
$3\lesssim \tilde {t}\lesssim 6$ corresponds to infiltration of the fine stratum. This uneven partitioning of flow is also reflected in the difference between the magnitudes of the superficial velocities 
$U_C=130\,\mathrm {\mu }{\rm m}\,{\rm s}^{-1}$ and 
$U_F=10\,\mathrm {\mu }{\rm m}\,{\rm s}^{-1}$ in the coarse and fine strata, respectively, corresponding to a ratio of 
$U_F/U_C=0.075$. We observe similar behaviour with our Newtonian control, which produces a similar ratio of 
$(U_F/U_C)_{0}=0.063$ even at a larger imposed flow rate 
$Q=35\,\mathrm {ml}\,\mathrm {h}^{-1}$ (movie 2). Hence, at low 
${Wi}_I$, polymer solutions recapitulate the uneven partitioning of flow across strata that is characteristic of Newtonian fluids.
Next, we repeat the same experiment as in figure 2(a i) at a larger flow rate of 
$Q=25\,\mathrm {ml}\,\mathrm {h}^{-1}$ – corresponding to a larger 
${Wi}_I=2.7$. Surprisingly, under these conditions, the partitioning of flow is markedly less uneven (figure 2(a ii), movie 3). These observations are reflected in the dye breakthrough curve, as well: the previously distinct jumps in the concentration 
$\tilde {C}$ begin to merge, as shown by comparing the blue and green lines in figure 2(b). Indeed, the ratio between the superficial velocities in the fine and coarse strata 
$U_F/U_C=0.16$, is 
${\sim }3{\times }$ larger than in the laminar baseline given by the Newtonian control and the low 
${Wi}_I=1.4$ solution tests. Therefore, to quantify this net improvement in flow homogenization, we normalize the velocity ratio by its Newtonian value, 
$\tilde {U}_F/\tilde {U}_C\equiv (U_F/U_C)/(U_F/U_C)_{0}=2.6$. This improvement in the flow homogenization is weaker at an even larger flow rate 
$Q=45\,\mathrm {ml}\,\mathrm {h}^{-1}$ (corresponding to 
${Wi}_I=3.3$), as shown in figure 2(a iii), the red line in figure 2(b), and in movie 4; the corresponding velocity ratio is 
$\tilde {U}_F/\tilde {U}_C=1.7$. Taken together, our observations demonstrate that polymer additives can help mitigate uneven partitioning of flow in a stratified porous medium – but that this effect is optimized at intermediate 
${Wi}_I$.
Why does this flow homogenization arise? As described in § 1, for the initially laminar flow, the coarser stratum experiences a higher interstitial flow speed – and therefore shear rate – as prescribed by the differential partitioning of flow across strata following Darcy's law. Thus, the locally defined Weissenberg number is larger in the coarse stratum, which we expect leads to an earlier onset of the elastic instability in this stratum. The corresponding increase in the resistance to flow through the coarse stratum would then redirect fluid to the fine stratum, helping to homogenize the flow. We test this expectation using our ‘continuous’ imaging protocol, which enables us to directly image the flow at the pore scale within the stratified microfluidic assembly. At the intermediate 
${Wi}_I=2.7$ – at which the flow homogenization is optimized – all pores observed in the fine stratum exhibit laminar flow that is steady over time (movie 5; representative pore shown in figure 2c i,d i). By contrast, 20 % of the pores observed in the coarse stratum exhibit strong spatial and temporal fluctuations in the flow (figure 2e). The fluid streamlines continually cross and vary over time, indicating the emergence of the elastic instability, as shown in movie 5 and in figure 2(c iii,d iii) for a representative pore – consistent with our expectation. These random streamline fluctuations are similar to those observed for this instability in a homogeneous medium (Browne & Datta Reference Browne and Datta2021). To highlight the regions of unstable flow, figure 2(c) also includes an overlay in red showing the standard deviation of the fluctuations in the fluorescent intensity over time; figure 2(d) shows the corresponding root mean square velocity fluctuation, which confirms that the flow is unstable in those same regions. At the even larger 
${Wi}_I=3.3$ – at which the improvement in flow homogenization is weaker – a larger fraction of pores in both strata exhibit unstable flow (movie 6; figure 2c ii,iv,d ii,iv, figure 2e). These results thus indicate that macroscopic flow homogenization is indeed linked to the onset of the elastic flow instability in the coarse stratum at sufficiently large 
${Wi}_I$, but is mitigated by the additional onset of the instability in the fine stratum at even larger 
${Wi}_I$.
3.2. Flow velocity fluctuations generated by the elastic flow instability lead to an increase in the apparent viscosity
To quantitatively understand the link between pore-scale differences in this flow instability and macro-scale differences in superficial velocity between strata, we consider the resistance to flow in the distinct strata at different 
${Wi}_I$. In particular, we model the strata as parallel fluidic ‘resistors’ – that is, we treat each stratum as a homogeneous porous medium (e.g. coarse 
$C$ or fine 
$F$), with the two hydraulically connected only at the inlet and outlet with fully developed flow in each. Because the time-averaged pressure drop 
$\langle \Delta P\rangle _t$ is equal across both strata, the imposed constant volumetric flow rate 
$Q$ must partition into the coarse and fine strata with flow rates 
$Q_C$ and 
$Q_F$, respectively, in proportion to their individual flow resistances via Darcy's law
\begin{equation} \frac{\langle\Delta P\rangle_t}{L}=\frac{\eta_{{app},C}Q_C}{k_CA_C} =\frac{\eta_{{app},F}Q_F}{k_FA_F}. \end{equation}
Following our previous study of this elastic flow instability in a homogeneous porous medium (Browne & Datta Reference Browne and Datta2021), we combine macro-scale pressure drop measurements with pore-scale flow visualization to determine and validate a model for the 
$\eta _{{app},i}$ of each stratum in isolation. We then use this model to deduce the apparent viscosity and uneven partitioning of flow within a stratified medium.
To do so, we measure the time-averaged pressure drop 
$\langle \Delta P\rangle _t$ at different volumetric flow rates 
$Q$ across each microfluidic assembly. We use Darcy's law to determine the corresponding 
$\eta _{{app}}$, which we plot as a function of 
${Wi}_I$ in figure 3(a); the data for the coarse medium are taken from Browne & Datta (Reference Browne and Datta2021). As expected, at small 
${Wi}_I\lesssim 2.6$, the apparent viscosity 
$\eta _{{app}}$ is given by the bulk solution shear viscosity 
$\eta (\dot {\gamma }_I)$, indicated by the red dashed line. However, above a threshold 
${Wi}_c=2.6$, 
$\eta _{{app}}$ rises sharply, paralleling previous reports (Marshall & Metzner Reference Marshall and Metzner1967; James & McLaren Reference James and McLaren1975; Durst et al. Reference Durst, Haas and Kaczmar1981; Durst & Haas Reference Durst and Haas1981; Kauser et al. Reference Kauser, Dos Santos, Delgado, Muller and Saez1999; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016). Both the homogeneous coarse (dark blue circles) and fine (light blue circles) media exhibit a similar dependence of 
$\eta _{{app}}$ on 
${Wi}_I$ – indicating that for our geometrically similar packings, 
$\eta _{{app}}({Wi}_I)$ does not depend on grain size 
$d_p$.
Figure 3. The elastic flow instability produces a similar increase in the macroscopic flow resistance for homogeneous porous media of different permeabilities. (a) Points show the apparent viscosity, normalized by the shear viscosity of the bulk solution, obtained using macroscopic pressure drop measurements. The apparent viscosity increases above a threshold 
${Wi}_I$ due to the onset of the elastic flow instability. Measurements for two different homogeneous media with distinct bead sizes and permeabilities (different colours) show similar behaviour. Grey line shows the predicted apparent viscosity using our power balance ((3.2), neglecting strain history) and the measured power-law fit to 
$\langle \chi \rangle _{t,V}$ shown in (b), with no fitting parameters; the uncertainty associated with the fit yields an uncertainty in this prediction, indicated by the shaded region. At the largest 
${Wi}_I$, the apparent viscosity eventually converges back to the shear viscosity, reflecting the increased relative influence of viscous dissipation from the base laminar flow. (b) Points show the rate of added viscous dissipation due to unstable flow fluctuations averaged over the medium, 
$\langle \chi \rangle _{t,V}$, measured from flow visualization. The dissipation sharply increases above the onset of the instability and is not sensitive to the bead size. Error bars represent one standard deviation between pores. We fit the data using an empirical power-law relationship 
${\sim }({Wi}_I/{Wi}_c-1)^{2.4}$ above the macroscopic threshold 
${Wi}_c=2.6$, shown by the grey line; the shaded region shows the error in the power-law fit.
To model this dependence of 
$\eta _{{app}}$ on 
${Wi}_I$, we directly image the pore-scale flow in each homogeneous microfluidic assembly with confocal microscopy using our ‘intermittent’ imaging protocol. We previously reported these measurements solely for the homogeneous coarse medium (Browne & Datta Reference Browne and Datta2021); thus, we first summarize these results. At small 
${Wi}_I<2.6$, the flow is laminar in all pores. Above the threshold 
${Wi}_c=2.6$, the flow in some pores becomes unstable, exhibiting strong spatio-temporal fluctuations. At progressively larger 
${Wi}_I$, an increasing fraction of the pores becomes unstable. To directly compute the added viscous dissipation arising from these flow fluctuations, we use our PIV measurements to determine the fluctuating component of the strain rate tensor 
$\boldsymbol{\mathsf{s}}'=(\boldsymbol {\nabla }\boldsymbol {u}'+\boldsymbol {\nabla }\boldsymbol {u}'^{\mathrm {T}})/2$. The rate of added viscous dissipation per unit volume arising from these flow fluctuations is then given directly by 
$\langle \chi \rangle _t=\eta \langle \boldsymbol{\mathsf{s}}':\boldsymbol{\mathsf{s}}'\rangle _t$, which can be estimated from the measured 
$x$–
$y$ velocity field (Sharp, Kim & Adrian Reference Sharp, Kim and Adrian2000; Delafosse et al. Reference Delafosse, Collignon, Crine and Toye2011). As anticipated, the overall rate of added dissipation per unit volume 
$\langle \chi \rangle _{t,V}$ determined by averaging 
$\langle \chi \rangle _t$ across all imaged pores increases with 
${Wi}_I$ above the threshold 
${Wi}_c=2.6$ (figure 3(b), dark blue circles) as a greater fraction of pores becomes unstable.
Next, we repeat this procedure in the homogeneous fine medium (figure 3(b), light blue circles). Intriguingly, the overall rate of added dissipation per unit volume 
$\langle \chi \rangle _{t,V}$ does not significantly vary between media. Additionally measuring 
$\langle \chi \rangle _{t,V}$ using our ‘continuous’ imaging protocol in the homogeneous fine medium further corroborates this agreement (figure 3(b), light blue squares). We speculate that this collapse reflects that flow fluctuations do not have a characteristic length scale (Browne & Datta Reference Browne and Datta2021); further studies of the influence of confinement on 
$\langle \chi \rangle _{t,V}$ will be a useful direction for future work. Our data indicate that, for the experiments reported here, differences in grain size between homogeneous porous media are well-captured by 
${Wi}_I$. We therefore fit all the data by the single empirical relationship 
$\langle \chi \rangle _{t,V}=A_x({Wi}_I/{Wi}_c-1)^{\alpha _x}$, with 
$A_x=176\pm 1\,\mathrm {W}\,\mathrm {m}^{-3}$, 
$\alpha _x=2.4\pm 0.3$, and 
${Wi}_c=2.6$, shown by the grey line in figure 3(b).
Finally, we follow our previous work (Browne & Datta Reference Browne and Datta2021) to quantitatively link the pore-scale flow fluctuations generated by the elastic flow instability to 
$\eta _{{app}}({Wi}_I)$. The power density balance for viscous-dominated flow relates the rate of work done by the fluid pressure 
$P$ to the rate of viscous energy dissipation per unit volume: 
$-\boldsymbol {\nabla }\boldsymbol {{\cdot }} P\boldsymbol {u}=\boldsymbol {\tau }:\boldsymbol {\nabla }\boldsymbol {u}$, where 
$\tau$ and 
$\boldsymbol {\nabla }\boldsymbol {u}$ are the stress and velocity gradient tensors, respectively. Averaging this equation over time 
$t$ and the entire volume 
$V$ of a given porous medium, and decomposing the velocity field into the sum of a base temporal mean and an additional component due to velocity fluctuations, then yields
\begin{equation}
\frac{\langle\Delta P\rangle_{t}}{\Delta
L}\equiv\frac{\eta_{{app}}
(Q/A)}{k}\approx{\underbrace{\frac{\eta(\dot{\gamma}_{I})
(Q/A)}{k}} _{\text{Darcy's law}}}\,
+{\underbrace{\frac{\langle\chi\rangle_{t,V}}{(Q/A)}}_{\text{Fluctuations}}}
+\left\{\begin{array}{@{}c@{}}\text{Strain}\\ \text{history}\\
\text{effects}\end{array}\right\}.
\end{equation}
The first term on the right-hand side of (3.2) represents Darcy's law for the base temporal mean of the flow. The second term reflects the added viscous dissipation by the solvent induced by the unstable flow fluctuations; our previous measurements (Browne & Datta Reference Browne and Datta2021) indicate that this dissipation does not vary appreciably along the imposed flow direction. The final term represents additional contributions arising from the full dependence of stress 
$\tau$ on polymer strain history in three dimensions (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987), which is currently inaccessible in our experiments. However, our previous measurements in the homogeneous course medium (Browne & Datta Reference Browne and Datta2021) indicate that this final term is relatively small for the range of 
${Wi}_I$ considered here, because the flow is quasi-steady and polymers do not accumulate appreciable Hencky strain over a duration of one polymer relaxation time 
$\lambda$. Therefore, for simplicity, we consider just the first two terms, which yields the grey line in figure 3(a); the shaded region indicates the uncertainty in this model arising from the empirical fit in figure 3(b). Our modelled 
$\eta _{{app}}({Wi}_I)$ thereby obtained from the pore-scale imaging shows excellent agreement with the 
$\eta _{{app}}$ obtained from the macro-scale pressure drop measurements (symbols) for both homogeneous media, without using any fitting parameters, for 
${Wi}_I\lesssim 4$. The slight discrepancies at larger 
${Wi}_I$ suggest that strain history effects play a non-negligible role in this regime. Nevertheless, as a first approximation, we use the 
$\eta _{{app}}({Wi}_I)$ modelled using (3.2) (neglecting the last term describing strain history) to deduce the apparent viscosity 
$\eta _{{app},i}$ within each stratum in (3.1).
3.3. Parallel-resistor model recapitulates experimental measurements of apparent viscosity and uneven flow partitioning
We next incorporate our model for the apparent viscosity 
$\eta _{{app},i}({Wi}_I)$ in the parallel-resistor model of a stratified medium described previously. Specifically, for a given imposed total flow rate 
$Q$, which corresponds to a given 
${Wi}_I$, we numerically solve (3.1) and (3.2) (neglecting the last term) along with mass conservation (
$Q=Q_F+Q_C$) to obtain the apparent viscosity 
$\eta _{{app}}({Wi}_I)$ for the entire stratified system.
To validate this approach, we first compute 
$\eta _{{app}}({Wi}_I)$ for the case of 
$\tilde {k}=9$ and 
$\tilde {A}=1$, which describes the stratified microfluidic assembly used in our experiments. Notably, the model shows a similar threshold 
${Wi}_c=2.6$ and overall shape of 
$\eta _{{app}}({Wi}_I)$ as in the homogeneous case, as shown by the blue line in figure 4(a) – suggesting that stratification does not appreciably alter the macroscopic flow resistance. Indeed, we find good agreement between this model prediction and our experimentally determined 
$\eta _{{app}}$, obtained from pressure drop measurements across the stratified microfluidic assembly, as shown by the blue circles in figure 4(a).
Figure 4. Parallel-resistor model captures the key features of experimentally measured apparent viscosity and uneven flow partitioning in stratified media. (a) Points show the normalized apparent viscosity measured for a stratified microfluidic assembly, indicating that it shows a similar increase above the onset of the elastic flow instability. Blue line shows the predicted apparent viscosity using our parallel-resistor model with no fitting parameters. Grey line shows the corresponding prediction for a homogeneous medium. Left and right arrows show 
${Wi}_I=2.7$ and 
$3.3$, at which only the coarse stratum or both strata are unstable in figure 2(c,d), respectively. The downward and upward triangles indicate the 
${Wi}_I$ at which each stratum becomes unstable. (b) Points show the ratio of superficial velocities in each stratum, normalized by the Newtonian value, measured for a stratified Hele-Shaw assembly; 
$\tilde {U}_{F}/\tilde {U}_C$ increases above the onset of the instability in the coarse stratum, indicating flow homogenization, and then decreases above the onset of the instability in the finer stratum as well, indicating that flow homogenization is mitigated. Teal line shows the prediction from our parallel-resistor model, which captures this non-monotonic behaviour.
This model also enables us to predict the onset of the elastic flow instability in the different strata at different values of the macroscopic 
${Wi}_I$. As demonstrated by the experiments on homogeneous media (figure 3), a given stratum becomes unstable when the local Weissenberg number exceeds the threshold 
${Wi}_c=2.6$. However, because of the difference in the permeabilities of the strata, flow partitions unevenly across them, causing different strata to reach this threshold at different imposed macroscopic 
${Wi}_I$. For small 
${Wi}_I$, the flow is slower in the fine stratum, with the ratio of superficial velocities given by the Newtonian value 
$(U_F/U_C)_0=0.075$. As a result, the model predicts that the coarse stratum becomes unstable at a smaller value of the macroscopic 
${Wi}_{c,C}=2.3$ (downward triangles in figure 4), and that the fine stratum becomes unstable at an even larger 
${Wi}_{c,F}=2.8$ (upward triangles). This prediction is in excellent agreement with our experimental pore-scale observations (figure 2c–e) that at 
${Wi}_{I}=2.7$ (left arrow in figure 4a), only the coarse stratum is unstable, while at a larger 
${Wi}_{I}=3.3$ (right arrow), both strata are unstable.
The model also reproduces and sheds light on the physics underlying the flow homogenization induced by the elastic flow instability, as we observed experimentally in the stratified Hele-Shaw assembly (figure 2a,b). For this case of 
$\tilde {k}=26$ and 
$\tilde {A}=1$, we use the model to compute the normalized ratio of superficial velocities 
$\tilde {U}_F/\tilde {U}_C\equiv (U_F/U_C)/(U_F/U_C)_0$ as a function of 
${Wi}_I$. The model prediction is shown by the line in figure 4(b). As expected, with increasing 
${Wi}_I$, the onset of the instability in the coarse stratum increases the resistance to flow in this stratum, redirecting fluid toward the fine stratum and thereby homogenizing the uneven flow across the entire medium – as indicated by the rapid increase in 
$\tilde {U}_F/\tilde {U}_C$ above 
${Wi}_{c,C}=2.3$ (downward triangle). However, this homogenization only arises in a window of flow rates: at even larger 
${Wi}_{I}>{Wi}_{I,F}=2.8$ (upward triangle), 
$\tilde {U}_F/\tilde {U}_C$ peaks and continually decreases, reflecting the onset of the instability in the fine stratum as well. While we do not expect perfect quantitative agreement with the experiments, given the assumptions and approximations made in our model, the experimental measurements show similar behaviour: as shown by the circles in figure 4(b), 
$\tilde {U}_F/\tilde {U}_C$ initially rises for 
${Wi}_I>{Wi}_{c,C}=2.3$, and then continues to decrease as 
${Wi}_I$ exceeds 
${Wi}_{c,F}=2.8$.
Thus, despite its simplicity, the parallel-resistor model of a stratified medium (3.1) that explicitly incorporates the increase in flow resistance generated by the elastic flow instability in each stratum (3.2) captures our key experimental findings: (i) the form of the macroscopic 
$\eta _{{app}}({Wi}_I)$ describing the entire medium, (ii) the differential onset of the instability in the different strata at varying 
${Wi}_I$ and (iii) the corresponding window of 
${Wi}_I$ within which the uneven flow across strata is homogenized. Having thereby validated the model, we next use it to further examine how the instability may homogenize fluid transport in stratified porous media having a broader range of permeability and area ratios, 
$\tilde {k}\equiv k_{C}/k_{F}$ and 
$\tilde {A}\equiv A_{C}/A_{F}$, respectively, than currently accessible in the experiments.
3.4. Geometry dependence of flow homogenization
How do the onset of and extent of homogenization imparted by this elastic flow instability depend on the geometric characteristics of a stratified porous medium? To address this question, we use our model to probe how the overall apparent viscosity 
$\eta _{{app}}({Wi}_I)$ and the flow velocity ratio 
$\tilde {U}_F/\tilde {U}_C({Wi}_I)$ depend on 
$\tilde {k}$ and 
$\tilde {A}$.
The measurements shown in figure 4 indicate that, despite the structural heterogeneity and uneven partitioning of the flow in a stratified medium, 
$\eta _{{app}}({Wi}_I)$ is not strongly sensitive to stratification; instead, it follows a similar trend to that of a homogeneous medium (
$\tilde {k}=1$). The model further supports this finding; with increasing 
$\tilde {k}$ (fixing 
$\tilde {A}=1$), the profile of 
$\eta _{{app}}({Wi}_I)$ shifts ever so slightly to smaller 
${Wi}_I$, eventually converging to the same final profile for 
$\tilde {k}\gg 100$, as shown in figure 5(a).
Figure 5. Geometry dependence of the apparent viscosity and uneven flow partitioning in a stratified medium, as predicted by our parallel-resistor model. (a,b) Different colours show the predictions of the parallel-resistor model for stratified media with varying ratios of the strata permeabilities, 
$\tilde {k}$, holding the area ratio fixed at 
$\tilde {A}=1$. The apparent viscosity (a) only shifts slightly to smaller 
${Wi}_{I}$ with increasing 
$\tilde {k}$, eventually converging for 
$\tilde {k}\gg 100$. The extent of flow homogenization generated by the elastic flow instability, quantified by the ratio of superficial velocities (b), does increase with increasing 
$\tilde {k}$. Optimal flow homogenization is indicated by the open circles at 
${Wi}_I={Wi}_I^{peak}$ with a velocity ratio 
$(\tilde {U}_F/\tilde {U}_C)^{peak}$. Inset to (a) shows the critical 
${Wi}_{I}$ at which each stratum becomes unstable; the window between the two values increases with increasing 
$\tilde {k}$. (c,d) Similar results to (a,b), but for stratified media with varying strata area ratios, 
$\tilde {A}$, holding the permeability ratio fixed at 
$\tilde {k}=9$. Inset to (c) shows the critical 
${Wi}_{I}$ at which each stratum becomes unstable; the window between the two values decreases with increasing 
$\tilde {A}$. Insets to (d) show the variation of the optimal 
${Wi}_I^{peak}$ and 
$(\tilde {U}_F/\tilde {U}_C)^{peak}$ with 
$\tilde {k}$, for different 
$\tilde {A}$. The data for different 
$\tilde {A}$ trivially collapse due to the definition of the superficial velocity.
However, the onset of the elastic flow instability in the different strata does vary with increasing 
$\tilde {k}$ (inset of figure 5a): 
${Wi}_{c,C}$ correspondingly shifts to slightly smaller 
${Wi}_I$, while 
${Wi}_{c,F}$ progressively shifts to larger 
${Wi}_I$, reflecting the increasingly uneven partitioning of the flow imparted by increasing permeability differences. As a result, the strength of the flow homogenization generated by the instability, as well as the window of 
${Wi}_I$ at which it occurs, increases with 
$\tilde {k}$ (figure 5b). This phenomenon is optimized at the peak position indicated by the open circles, which occur at 
${Wi}_I={Wi}_I^{peak}$ with a flow velocity ratio 
$(\tilde {U}_F/\tilde {U}_C)^{peak}$. We therefore summarize our results by plotting both quantities as a function of 
$\tilde {k}$ (dark blue lines, insets to figure 5d). Again, both increase until 
$\tilde {k}\approx 400$. For even larger 
$\tilde {k}$, 
${Wi}_I^{peak}$ plateaus at 
$\approx 3.7$, while 
$(\tilde {U}_F/\tilde {U}_C)^{peak}$ plateaus at 
$\approx 5.5$. This behaviour reflects the non-monotonic nature of our model for 
$\eta _{{app},i}({Wi}_I)$; at such large permeability ratios, the coarse stratum reaches its maximal value of 
$\eta _{{app},C}$ at 
${Wi}_I<{Wi}_{c,F}$, maximizing the extent of flow redirection to the fine stratum generated by the instability in the coarse stratum. These physics are also reflected in the values of 
${Wi}_I^{peak}$ and 
${Wi}_{c,F}$ (open circles and filled upward triangles in figure 5b, respectively); while the two match for small 
$\tilde {k}$, 
${Wi}_I^{peak}$ becomes noticeably smaller than 
${Wi}_{c,F}$ for 
$\tilde {k}\gtrsim 400$.
Similar results arise with varying 
$\tilde {A}$ (fixing 
$\tilde {k}=9$), as shown in figure 5(c,d). Here, 
$\tilde {A}<1$ and 
$\tilde {A}>1$ describe the case in which a greater fraction of the medium cross-section is occupied by the fine or coarse stratum, respectively; the limits of 
$\tilde {A}\rightarrow 0$ and 
$\rightarrow \infty$ therefore represent a non-stratified homogeneous medium. While stratification again does not strongly alter 
$\eta _{{app}}({Wi}_I)$, we find that 
${Wi}_{c,C}$, 
${Wi}_{c,F}$, and 
${Wi}_I^{peak}$ increase with 
$\tilde {A}$. Furthermore, 
$(\tilde {U}_F/\tilde {U}_C)^{peak}$ does not depend on 
$\tilde {A}$, since the superficial velocity incorporates cross-sectional area by definition. Taken together, these results provide quantitative guidelines by which the macroscopic flow resistance, as well as the onset and extent of flow homogenization, can be predicted for a porous medium with two parallel strata of a given geometry.
3.5. Extending the model to porous media with even more strata
As a final demonstration of the utility of our approach, we extend it to the case of a porous medium with 
$n$ parallel strata, each indexed by 
$i$. To do so, we again maintain the same pressure drop across all the different strata (3.1), with the apparent viscosity 
$\eta _{{app},i}$ in each given by (3.2), and numerically solve these 
$n-1$ equations constrained by mass conservation, 
$Q=\varSigma _{i=1}^{n} Q_i$.
As an illustrative example, we consider 
$n=5$ with the different stratum permeabilities chosen from a log-normal distribution, as is often the case in natural settings (Freeze Reference Freeze1975): 
$k_i\in \{79,51,36,26,17\}\,\mathrm {\mu }\mathrm {m}^2$. To characterize the flow redirection between strata at varying overall 
${Wi}_I$, we focus on the ratio of the superficial velocity 
$U_i$ in each stratum and the macroscopic superficial velocity 
$U\equiv Q/A$, normalized by the value of this ratio for a Newtonian fluid: 
$\tilde {U}_{i}/\tilde {U}\equiv (U_{i}/U)/(U_{i}/U)_0$. Hence, larger (smaller) values of 
$\tilde {U}_{i}/\tilde {U}$ indicate that fluid is being redirected to (from) a given stratum 
$i$. Consistent with our previous results, the coarsest stratum becomes unstable at the smallest 
${Wi}_I$ (dark purple line in figure 6a), redirecting fluid to the other strata – as indicated by the reduction in 
$\tilde {U}_{i}/\tilde {U}$ for 
$k_i=79\,\mathrm {\mu }\mathrm {m}^2$ as 
${Wi}_I$ increases above 
$\approx 2.4$, and the concomitant increase in 
$\tilde {U}_{i}/\tilde {U}$ for the other strata (blue to light green lines). Each progressively finer stratum then becomes unstable at progressively larger 
${Wi}_I$, as indicated by the upward triangles, redirecting fluid from it to the other strata. Thus, as with the case of 
$n=2$ examined previously, the flow homogenization generated by the elastic flow instability arises only in a window of 
${Wi}_I$.
Figure 6. Model predictions for a porous medium with five distinct strata. (a) Different colours show the predicted ratio between the superficial velocity in each stratum and the macroscopic superficial velocity, normalized by the value of this ratio for a Newtonian polymer-free solvent. The coarsest stratum (dark purple) becomes unstable at the smallest 
${Wi}_I$, following by the next coarsest (dark blue) and so on – causing flow to be redirected to the finer strata and the uneven flow across different strata to be homogenized. At even larger 
${Wi}_I$, all the strata become unstable and the resulting flow homogenization is mitigated. (b) Predicted breakthrough curves for the polymer solution at 
${Wi}_I=3.2$ (light green) as well as the Newtonian polymer-free solvent at the same flow rate (dark green). At this intermediate 
${Wi}_I$, the elastic flow instability homogenizes the uneven flow across strata; as a result, rapid breakthrough in the coarsest strata is slowed (left arrow), and slow breakthrough in the finest strata is hastened (right arrow), smoothing the overall breakthrough curve. Inset shows the macroscopic effective longitudinal dispersivity, normalized by its value for the Newtonian polymer-free solvent at the same volumetric flow rate. Here, 
$K_l$ and 
$K_{l,0}$ differ slightly at low 
${Wi}_I$ because of the modest shear thinning in the polymer solution, which increases the uneven partitioning of flow uniformly before the onset of unstable flow. For a window of 
$2.4\lesssim {Wi}_I\gtrsim 4.5$, the normalized dispersivity is smaller than one, indicating more uniform scalar transport due to the homogenized flow resulting from the instability.
As a final illustration of this point, we compute the corresponding breakthrough curve of a passive scalar, 
$\tilde {C}(t)$, given that such curves are commonly used to characterize transport in porous media for a broad range of applications. To do so, for a given stratum 
$i$ with 
$U_i$ determined from our parallel-resistor model, we use the foundational model of Perkins & Johnston (Reference Perkins and Johnston1963) as an example to compute
\begin{equation} C_i(t)=0.5\left[1-\mathrm{erf}\left(\frac{1-t/t_{PV}}{2\sqrt{K_{l,i}/U_iL}\sqrt{t/t_{PV}}}\right)\right]. \end{equation}
This expression explicitly incorporates the dispersion of a passive scalar being advected by the flow via the longitudinal dispersivity 
$K_{l,i}$, which depends on the scalar diffusion coefficient 
$D$, the stratum tortuosity 
$\tau$ and the Péclet number characterizing scalar transport in a pore 
${Pe}=U_i d_{p,i}/D$; in particular, 
$K_{l,i}=D(1/\tau +0.5{Pe}^{1.2})$ when 
${Pe}<605$ and 
$K_{l,i}=D(1/\tau +1.8{Pe})$ when 
${Pe}>605$ (Woods Reference Woods2015). The overall breakthrough curve is then given by 
$C(t)=\sum _{i}^{n}C_i(t)A_i/A$, which we normalize by its maximal value at 
$t\rightarrow \infty$ to obtain 
$\tilde {C}(t)$. For this illustrative example, we use values characteristic of small molecule solutes in natural porous media: 
$D=10^{-6}\,\mathrm {cm}^2\,\mathrm {s}^{-1}$, 
$\tau =2$ (Datta et al. Reference Datta, Chiang, Ramakrishnan and Weitz2013), and estimate 
$d_{p,i}$ from the stratum permeability using the Kozeny–Carman relation (Philipse & Pathmamanoharan Reference Philipse and Pathmamanoharan1993).
The resulting breakthrough curves 
$\tilde {C}(t)$ are shown in figure 6(b) for a fixed flow rate, chosen such that 
${Wi}_I=3.2$ for our polymer solution – just above the onset of the elastic flow instability in the finest stratum, at which we expect flow homogenization to be nearly optimized (figure 6a). For the case of the polymer-free Newtonian solvent, the flow partitions unevenly across the strata, leading to highly heterogeneous scalar breakthrough. As shown by the dark green line, coarser strata are infiltrated rapidly, leading to the rise in 
$\tilde {C}(t)$ at 
$t/t_{PV}\approx 0.4$. However, the considerably smaller flow speeds in the bypassed finer strata give rise to far slower breakthrough, leading to the subsequent jumps in 
$\tilde {C}(t)$ at longer times; as a result, 90 % of scalar breakthrough only occurs after 
$t/t_{PV}=2.5$ has elapsed. The polymer solution exhibits strikingly different behaviour: the breakthrough curve shown by the light green line is noticeably smoother, reflecting the flow homogenization imparted by the elastic flow instability. In this case, unstable flow hinders rapid infiltration in the coarser strata (right-pointing arrow at 
$t/t_{PV}\approx 0.6$), instead redirecting fluid to the finer strata (left-pointing arrow at 
$t/t_{PV}\approx 2$); as a result, 90 % of scalar breakthrough occurs 
$\approx 1.4\times$ faster, at 
$t/t_{PV}=1.8$.
This improvement in scalar breakthrough can also be described using an effective, macroscopic, stratum-homogenized longitudinal dispersivity 
$K_{l}$. Despite the complex shapes of breakthrough curves that commonly arise for stratified porous media due to uneven flow partitioning (e.g. dark green line in figure 6b), a standard practice is to fit the entire breakthrough curve to a single error function (Lake & Hirasaki Reference Lake and Hirasaki1981) and thereby extract 
$K_{l}$. The dispersivity thereby determined from our computed breakthrough curves is shown in the inset to figure 6(b) for a broad range of 
${Wi}_I$. At small 
${Wi}_I$, 
$K_{l}$ matches that of a polymer-free Newtonian solvent 
$K_{l,0}$ at the same volumetric flow rate. Above 
${Wi}_I\approx 2.4$, at which the coarsest stratum becomes unstable, 
$K_{l}$ drops relative to the Newtonian value – indicating more uniform scalar breakthrough due to flow homogenization. The effective dispersivity continues to decrease as an increasing number of strata become unstable, further homogenizing the flow and causing scalar breakthrough to become more uniform. The effective dispersity is ultimately minimized at the optimal 
${Wi}_I\approx 3.2$. Increasing 
${Wi}_I$ further causes 
$K_{l}/K_{l,0}$ to then increase, eventually reaching 1 at 
${Wi}_I\approx 4.5$ – again reflecting the fact that the flow homogenization generated by the elastic flow instability arises in the window of 
$2.4\lesssim {Wi}_I\lesssim 4.5$.
4. Conclusions
The work described here provides the first, to our knowledge, characterization of a purely elastic flow instability in stratified porous media. Our experiments combining flow visualization with pressure drop measurements revealed that the instability arises at different flow rates, corresponding to different 
${Wi}_I$, in different strata. Uneven partitioning of flow into the higher-permeability strata causes them to become unstable at smaller 
${Wi}_I$ – redirecting the flow towards the lower-permeability strata, thereby helping to homogenize the flow across the entire medium. At even larger 
${Wi}_I$, the lower-permeability strata become unstable as well, suppressing this flow redirection – leading to a window of flow rates at which this homogenization arises.
We elucidated the physics underlying this behaviour using a minimal parallel-resistor model of a stratified medium that explicitly incorporates the increase in flow resistance generated by the elastic flow instability in each stratum. Despite the simplicity of the model, it captures the macroscopic resistance to flow through the entire medium, the differential onset of the instability in the different strata at varying 
${Wi}_I$, and the corresponding window of 
${Wi}_I$ within which the uneven flow across strata is homogenized, as found in the experiments. Taken together, our work thus establishes a new approach to homogenizing fluid and passive scalar transport in stratified porous media – a critical requirement in many environmental, industrial and energy processes.
This study focused on a single polymer solution formulation as an illustrative example. However, the threshold 
${Wi}_c$ at which the instability arises, and the corresponding excess flow resistance 
$\langle \chi \rangle _{t,V}$, likely depend on the solution rheology (through e.g. polymer concentration, molecular weight and solvent composition). The relative importance of the full polymer strain history in three dimensions, neglected here for simplicity, may also play a non-negligible role for different formulations and at large 
${Wi}_I$; indeed, while we use the specific functional form of 
$\eta _{{app}}$ given by (3.2), it is unclear how far this model can be extrapolated past 
${Wi}_I\gtrsim 4$. Incorporating these additional complexities into our analysis will be an important direction for future work.
Nevertheless, the theoretical framework established here provides a way to develop quantitative guidelines for the design of polymeric solutions and fluid injection strategies, given a stratified porous medium of a particular geometry. We therefore anticipate it will find use in diverse applications – particularly those that seek to balance the competing demands of minimizing the macroscopic resistance to flow (quantified by 
$\eta _{{app}}$) and maximizing flow homogenization (quantified by 
$\tilde {U}_i$). Indeed, accomplishing this balance is a critical challenge in subsurface processes such as pump-and-treat remediation of groundwater, in situ remediation of groundwater aquifers using injected chemical agents, enhanced oil recovery and maximizing fluid–solid contact for heat transfer in geothermal energy extraction – for which uneven flow across strata is highly undesirable. Moreover, similar flows also play key roles in determining separation performance in filtration and chromatography, and improving heat and mass transfer in microfluidic devices. Thus, by deepening fundamental understanding of how the elastic flow instability can be harnessed to homogenize flow in stratified media, we expect our results to inform a broader range of applications.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.337.
Acknowledgements
It is our pleasure to acknowledge the Stone Lab for use of the rheometer.
Funding
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program (to C.A.B.) under grant no. DGE-1656466, a Camille Dreyfus Teacher-Scholar Award from the Camille and Henry Dreyfus Foundation (to S.S.D.), and partial support from Princeton University's Materials Research Science and Engineering Center under NSF grant no. DMR-2011750. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. C.A.B. was also supported in part by a Mary and Randall Hack Graduate Award from the High Meadows Environmental Institute and a Wallace Memorial Honorific Fellowship from the Graduate School of Princeton University.
Declaration of interests
The authors report no conflict of interest.
Author contributions
C.A.B. and S.S.D. designed the experiments; C.A.B. performed all experiments; C.A.B. and S.S.D. designed the theoretical model; C.A.B. and R.B.H. performed all theoretical analysis and numerical simulations of 2-strata media; C.A.B. and C.W.Z. performed all theoretical analysis and numerical simulations of 
$n$-strata media; C.A.B. and S.S.D. analysed all data, discussed the results and implications and wrote the manuscript; S.S.D. designed and supervised the overall project.
Data availability statement
All data are available in the manuscript and supplementary movies.
Open access
