MODELING FLOW IN VUGGY MEDIA Advisor Todd Arbogast

MODELING FLOW IN VUGGY MEDIA Advisor Todd Arbogast - Mathematics, University of Texas at Austin Diana Battefeld - Physics, Brown University Cheng Chen - Engineering, Northwestern University Abdullah Cihan - Engineering, University of Tennessee, Knoxville Daniel Pringle - Geophysics, University of Alaska, Fairbanks Masa Prodanovic – Comp. Applied Math, University of Texas at Austin Danail Vassilev - Computational Mathematics, University of Pittsburgh Hua Zhang - Agronomy, Louisiana State University Summer School In Geophysical Porous Media July 17 -28, 2006 Purdue University, West Lafayette, Indiana

Motivation ØEnvironmental: Understanding movement of contaminants in vuggy subsurface systems, pumping and protection of water aquifiers. ØOil Industry: Oil recovery - Help oil industries to manage oil reserves more efficiently. ØGeological: System unlike those studied normally by geologists. ØMathematical: Model fluid flow and transport in vuggy (dual porosity) media.

What is a Vug? Vug: an orifice in a porous medium relatively large compared to pore size Internal structure of Vug system: X-ray Computed Tomography (CT) Scans of Pipe Creek Reef sample

Where? This sample: Pipe Creek, Texas River basin Area Limestone matrix with caprinid fossils 4 -6 cm diameter Cretaceous caprinid fossil

Goal: To characterize the vuggy media’s pore space and to understand flow and tracer transport. We use CT scans, permeability measurements, and experimental tracer curves to achieve this goal. We will proceed as follows: I. CT Image analysis used to characterize pore space & connectivity II. Permeability computation from simplified Darcy flow models III. Tracer analysis what we can infer from the experiments

I CT IMAGE ANALYSIS

3 DMA-Rock software http: //www. ams. sunysb. edu/~lindquis/3 dma_rock/3 dma_rock. html Segmentation (Indicator Kriging, Oh and Lindquist ’ 99) Medial Axis (MA) extraction (Lee-Kashyap-Chu ’ 94) + some trimming Throat finding + pore partitioning Interphase areas, fluid blob characterization, pore scale saturation… Pore/throat characterization (distributions)

Medial Axis § Centrally located skeleton Ø Preserves topology and geometry of a 3 D object Ø a powerful search tool § Obtained by morphological thinning algorithms (3 DMA-Rock implements Lee. Kashyap-Chu '94. algorithm) § Branch cluster (node) – digitized version of a point (graph vertex) where medial axis paths meet § Medial axis path (link) – digitized version of a curve (graph link) § First applied to (simulated) porous media by Thovert/Sales/Adler ’ 93. and to sandstones by Spanne/Thovert/Jacquin/Lindquist/Jones/Adle r ’ 94. Terminology : Medial axis = skeleton= percolating backbone = deformation retract MA computation = morphological skeletonization = morphological thinning/erosion Grass-fire algorithm = burn algorithm = erosion algorithm = distance labeling

Objective: Pore-throat network Most throat finding algorithms minimize area A (and not hydraulic radius R)! R = A/P = (a. L 2)/(b. L)=(a/b)L = c. L (a, b, c constants reflective of shape) Throat: R 1 > R 2 < R 3 § Hydraulic radius of fluid R 1 = c 1 L 1 & R 2 = c 2 L 2, assume c 1=c 2 pathway cross-section: for close cross-sections, R = Area/Perimeter Then § Throat: a pathway cross R 1 < R 2 iff L 1 < L 2 iff section of minimal R § Pore: pore space opening c 1 L 12 < c 2 L 22 iff A 1 < A 2 enclosed by grain & throat surfaces

Throat Finding Approaches § Delaunay tesselation of sphere packs Ø “too many” pores obtained; merging close clusters needed Ø Finney '70, Mason '71, Bryant ’ 93. , Willson ’ 00. § Morphological thinning based pore partitioning Ø Baldwin/Sederman/Mantle/Alexander/Gladden '96. Ø widely used in soil sci. ; pore space partitioned, no throats identified § Maximal ball algorithm Ø Silin/Patzek 2003. , results in “ball-n-stick diagram” § Throat finding algorithms 1. Multi-orientation scanning (Zhao/Mac. Donald/Kwiecien '94. ) 2. Dijkstra based throats (Venkatarangan/Lindquist '99. ) and extensions (Prodanovic/Lindquist/Seright ‘ 05) MA based 3. Planar throats (Liang/Ioannidis/Chatzis '99. ) MA based 4. Wedge based throats (Shin/Lindquist '02. ) MA based

Sample Analysis

Does upscaling by coarsening work? § cut central 300 x 100 subset of the segmented image § original z-direction slices multiplied 3 times to match resolution in x- & ydirections 3003 § the image coarsened by using majority rule in the boxes of appropriate size 3003 Porosity 17. 8% Connectivity 83. 3% 1003 Porosity 17. 2% Connectivity 83. 0% 503 Porosity 16. 9% Connectivity 84. 7% Connectivity = percentage of pore space in the largest connected component

Vug Space Medial Axis (MA) 3003 § 1003 503 Small (MA) structures are destroyed by upscaling (above). Main fluid pathways are preserved: the shortest z-direction pathways are shown below. 503 opens up one more path!

Medial Axis Path Statistics

Sample Heterogeneity § Shortest paths across x-, y- and z- directions for 1003 sample § Average geometric tortuosity (ratio of actual length and side-side distance) for the shortest paths is ~1. 8 for all directions, and ~2. 0 for all possible paths across the sample

Vug-throat network § 3 DMA-Rock throat finding algorithms used to identify individual vugs vug surface MA path Branch cluster A vug from 3003 sample, 9983 mm 3 Relative positions of 10 largest vugs from 3003 sample

Vug-throat Network Statistics Small structures disappear in the upscaled image • If small structures are important for the flow, the simulation in upscaled images will not produce good results • Dead-end vugs (not included the above volume stats) can be identified and occupy 4% of the connected pore space component

II PERMEABILITY COMPUTATION

Flow and Transport Modeling Provide link between structure and experiment. Can system be understood in terms of an effective permeability? What controls tracer transport? Matrix (no vugs) k = 10 m. Darcy Sub-sample (~10 cm)3 k = 100 Darcy Flow Advection Diffusion Equation

Flow and Transport: Approaches 1. Medial axis and statistics 2. CT output structures Simple pipe flow ‘Darcy code’: modeling Simple pipe networks 2 D Lattice Boltzmann Medial-axis pipe networks 3. Dual porosity model fitting

Insights from Poiseuille Flow 1 Darcy (D) ~ 10 -12 m 2 Parallel Series Large k dominate Small k dominate

Insights from Poiseuille Flow 1 Darcy (D) ~ 10 -12 m 2 d = 10 cm kmx =10 m. D, kpipe = 45 D n=2: keff = 90 Darcy r = 1 mm Small constrictions can strongly control keff

Pipe Network Model Dead end v=0 P=1 P=0 INPUT: Medial axis path network NB: more usual to use pore/throat network, but we had medial axis data available first!

Pipe Network Formulation q 4 Poiseuille flow q 2 P 4 P 6 q 1 Mass balance P 1=1 q 6 P 7 P 2 q 3 P 5 q 8 P 8=0 q 7 q 5 Conjugate gradient method on normal equations: ATAx=ATb

Simulated Permeability Kzu Kxu Kzd Kxub Kxd Kyu Kyd Kzu Kzd Darcy 50× 50 a 1 x 106 2 x 106 8 x 106 1 x 106 100× 100 10500 30200 152 439 9630 839 300× 300 1320 1110 985 1400 3080 606 a Image resolutions. Medial axis flow paths are used for simulation. b Flow directions.

Simulated Permeability Kzu Kxu Kzd Kxub Kxd Kyu Kyd Kzu Kzd Darcy 50× 50 a 1 x 106 2 x 106 8 x 106 1 x 106 100× 100 10500 30200 152 439 9630 839 300× 300 1320 1110 985 1400 3080 606 a Image resolutions. Medial axis flow paths are used for simulation. b Flow directions. Upscaling dependence and Variability 2 -scale problem ! Voxel-sized restrictions to flow Vug scale ~ sample size

Darcy Code: Parssim Parallel subsurface simulator www. ices. utexas. edu/~arbogast/parssim/ IN: 3 D geometry + physical parameters OUT: (steady-state) velocity field, concentration field 48 x 48 voxels (i) Full sample, Coarse resolution L= 16. 4 cm) ( k = 416 Darcy (ii) Sub-samples, max resolution ( L = 2. 6 cm) k(s 1) = 1282 Darcy k(s 2) = 170 Darcy k(vug) = 107 D, k(mx) = 0. 01 D, ΔP = 104 Pa (0. 1 atm),

2 -D Lattice Boltzmann Modeling Flow and passive tracer transport Limited input file size: 100 x 100 voxels Q: How to project ‘very’ 3 D structure to 2 D ? A: Examine 1 -2 ‘connecting’ paths instead. Main flow Dead ends 2 D toy model made by tracing and projecting 2 paths

Flow Field

t = 0 (sec)

t = 85. 3 (sec)

t = 341. 2 (sec)

t = 1108. 9 (sec)

t = 1706 (sec)

t = 3412 (sec)

t = 6397. 5 (sec)

LBM Breakthrough Curve Reynolds number, Re = 0. 10 Inlet Permeability, k = 1. 1 x 106 Darcy Outlet Breakthrough Curves

Numerical Permeability Results Method Network Model (300)3 (Medial Axis) k [Darcy] 1400 ± 800 PARSSIM upscaled subsample 1 subsample 2 416 1282 170 2 D Lattice Boltzmann 1 x 106 (? ) Within 1 -2 orders of magnitude of experiments: this was our target! 2 scale problem: really need Lsample > vugs and voxel < mm To finish: 3 D Lattice Boltzmann, pore-throat network system Darcy code: systematic examination of upscaling and k values. New numerical scheme in Arbogast and Brunson, submitted 2006.

III TRACER TRANSPORT

WHAT WE HAVE SEEN SO FAR: A. VUGS ARE INTERCONNECTED B. WE CAN EXTRACT AN EFFECTIVE PERMEABILITY C. APPROXIMATION: DARCY FLOW (only valid at larger scales) D. BETTER: POISEUILLE FLOW NOW LET’S SEE WHAT WE EXPECT FROM EXPERIMENTS

TOY MODEL: 3 D cellular vug-network (periodic)

HEURISTIC ARGUMENTS USING STOKES FLOW: causes for marker behaviour • Early Breakthrough: Tracer finds a pathway through the network • Multiple Plateaus: Tracer finds other (longer) ways through the network – result: new plateaus whenever a new way opens. • Abrupt Drop: Time drop correlated to end of injection. A time delay is due to cells on the main paths, acting as tracer reservoirs for a small amount of time. • Single Plateau: Weakly interconnected vugs are drained slowly. The large main volumes of vugs act as the reservoir for a longer time, leading to a plateau of constant height until they are drained. • Long Tail: Remaining tracer in thin arms is slowly washed out after main volumes are drained – diffusive process (exponential decay). Implement model on computer: test arguments!

TEST HYPOTHESIS (Experimentally/Numerically, 3 D) • Multiple plateaus: Close one of the main pathways and see if the early plateaus change. • Abrupt drop: Increase injection time – drop should be shifted by same time (the delay does not change) • Single Plateau: Inject the fluid for smaller time interval – no reservoir build up: Plateau should decrease in height. • Long tail: increase the diameter of the smallest arms in the toy model: exp. decay should be faster.

TEST HYPOTHESIS (Experimentally or Numerically, 3 D) • Multiple plateaus: experimentally- close one of the main pathways and see if the early plateaus change. From CT scan we see multiple pathways In Poiseuille flow model, plateaus are non-existent (no reservoirs).

Macroscopic Models Dual Porosity Model Mobile Immobile Dispersion Advection Diffusion between mobile and immobile region

Dual Porosity Model CXTFIT Diffusion between mobile and immobile regions important in transport. ( i. e. from vugs in/out of porous matrix ). Tracer data CXTFIT free fit

Stream Tube Model L 1 r 4 L 3 r 2 r 3 Li effective flow length ri effective radius Correlated random variables with known distributions Solute transport in each tube follows ADE Hydrodynamic dispersion Di vi - Poiseuille flow Lo: Length of the column Tortuosity ti=(L 0/Li)2 Longitudinal dispersivity a=0. 01 m Flux averaged effluent concentration

Parallel 4 -Tube Model Tracer data 4 tube model Generally good fit! Fit could be optimized. . No connectivity in model!

3 D models (outlook) • Easy: use periodic cell model as depicted, insert fluid, solve flow numerically (maybe analytic estimate – use cylinders, Pois. flow. ) • Better: Construct fractal Vug-network: bifurcate progressively to smaller branches (thereafter numerics as above). • Best: Add Monte Carlo to choose (around mean): - Length of arms before bifurcation - Angle of outgoing arms - Diameter after bifurcation of arms Goal: find effective permeability, check heuristic arguments, … Doable in 2 -weeks time: Single Parallelpiped Uniform/Non-Uniform Darcy. Flow 3 D Model.

References • T. Arbogast and H. L Lehr, Homogenization of a Darcy-Stokes system modeling vuggy porous media. , 2003. (Preprint available www. ices. utexas. edu/reports/2002. html, TICAM report 02 -44) • H. L. Lehr, Analysis of a Darcy-Stokes system modeling flow through vuggy porous media, Ph. D Thesis 2004 • W. B. Lindquist, 3 DMA-Rock, A Software Package for Automated Analysis of Rock Pore Structure in 3 -D computed Microtomography Images http: //www. ams. sunysb. edu/~lindquis/3 dma_rock/3 dma_rock. html • S. Succi, Lattice Boltzmann Equations for Fluid Dynamics and Beyond (Numerical Mathematics and Scientific Computation), 2001. • M. C. Sukop, D. T. Jr. Thorne, Lattice Boltzmann Modeling: An introduction for Geoscientists and Engineering, 2005 • D. Zhang and Q. Kang, Pore Scale simulation of solute transport in fractured porous media. , Geophysical Research Letters, 31, 2004