Density driven flow in porous media: How accurate are our models? Wolfgang Kinzelbach Institute for Hydromechanics and Water Resources Engineering Swiss Federal Institute of Technology,
Download ReportTranscript Density driven flow in porous media: How accurate are our models? Wolfgang Kinzelbach Institute for Hydromechanics and Water Resources Engineering Swiss Federal Institute of Technology,
Density driven flow in porous media: How accurate are our models? Wolfgang Kinzelbach Institute for Hydromechanics and Water Resources Engineering Swiss Federal Institute of Technology, Zurich, Switzerland Contents • Examples of density driven flow in aquifers • Equations – Formulation – Special features of density driven flows • Benchmarks – Analytical and exact solutions – Experimental benchmark: Grid convergence – Experimental benchmark: Fingering problem • Upscaling issues • Conclusions Density driven flows in groundwater resources management • Sea water intrusion • Salt water upconing under freshwater lenses (both on islands and inland) • Salt water fingering under playa lakes and saltpans • Flow around salt domes (nuclear waste repositories) • Brine injection • Leachate from waste deposits • Even the ordinary tracer experiment... Saltwater Intrusion Fresh water Salt water Formation of toe Fresh water Salt water Saltwater Upconing Fresh water 1 40 Salt water Example: Salt Water Upconing on Wei Zhou Island Thesis Li Guomin Freshwater Lens Upconing Alternative Extraction Strategies Salinity backflow from Chotts in Tunisia Oasis Salt fingers on islands in the Okavango Delta 200 km Islands and salt crusts in Delta Schematic cross section of an island Transpiration Trona saltcrust Evaporation gravity vs. upward flow Instability on the Islands WoodingNumber k f uET kf= 10-5 m/s, uET=10-8 m/s Critical Wooding Number: TEM Evaluation 20 Island-Transect 18 Critical Wooding Number 16 instable 14 12 10 8 6 stable 4 2 0 0 0.2 0.4 0.6 0.8 1 1.2 k in units of u/D m -1 1.4 1.6 1.8 2 TDS center (mg/l) TDS layer 2 norm. density Wooding (mg/l) contrast (-) No.(-) Thata-1 22425.00 1352.94 0.0147 14.74 Thata-3 Mosupatsela-1 9425.00 8645.00 878.59 1495.31 0.0060 0.0050 5.98 5.00 Mosupatsela-3 5200.00 1803.28 0.0024 2.37 Tshwene-1 Tshwene-2 19110.00 19110.00 1135.50 1135.50 0.0126 0.0126 12.57 12.57 Lebolobolo-1 Lebolobolo-2 13000.00 26000.00 1014.03 1887.59 0.0084 0.0169 8.38 16.86 Lebolobolo-3 Monyopi-1 7800.00 2665.00 886.25 1118.03 0.0048 0.0011 4.84 1.08 Monyopi-2 2665.00 747.22 0.0013 1.34 Atoll-3 Kwena-1 15600.00 9750.00 2172.15 1480.68 0.0094 0.0058 9.39 5.78 Montsana-1 1625.00 367.65 0.0009 0.88 Simulation of fingering 20 40 60 80 100 20 t=900 t=6000 t=8500 t=16800 t=32500 dddd t=66000 t=46500 t=2900 t=12400 t=25000 40 60 80 100 100 cmax=350 =235 =11 mg/l =54 =75 ccmax =350 mg/l =110 =30 max Flow in the vicinity of a salt dome Recharge Discharge Salt water fresh water interface Top of salt dome No density difference With density difference Basic Equations expressed in mass fraction c and pressure p • • • Darcy law • Dispersion tensor • • (n) (v) (c q )s vol t Mass balance total mass Mass balance salt (nc) ( cv n( D Dm )c) cq (cq ) svol t K v (p g) Dij Constitutive relationships Boundary conditions (many combinations) 1 v a vv ijkl k l k ,l (c, p) (c, p) e.g. c 1 c 1 1 cmax f cmax s 1 Possible simplification: Boussinesq approximation Features of density driven flow -Non-linearity -Consistency problem of boundary conditions -Rotational flow with closed streamlines -Plus all difficulties known from advectivedispersive transport Flow in porous media and rotation Darcy-flow in heterogeneous porous media is rotational Example: v kf But we still have: v k f h v 0 (v / k f ) 0 In density flow, rotation is non-trivial: closed streamlines k v (p g) k k For constant k/ v ( (p) ( ) g ) ( ) g Rotational when not parallel to g Numerical solution and testing of codes • • • • • Analytical solutions Exact solutions Inter-code comparison Experimental benchmarks Grid convergence All computations are made with d3f, a density flow model using unstructured grids, finite volume discretization, multigrid solver, error estimator, automatic local refinement/coarsening, parallel computing Idea of „exact“ solution (steady state) Pressure equation L1 ( p( x, z), c( x, z) 0 Salt mass fraction equation L2 (c( x, z), p( x, z)) 0 Assume any differentiable functions p(x,z), c(x,z) Assume any domain Assign function values as first kind boundary conditions on boundary of that domain Plug functions into flow equations Pressure L1 ( p( x, z), c( x, z)) f ( x, z) Salt mass fraction L2 (c( x, z), p( x, z)) g ( x, z) Right-hand sides are not zero: They are interpreted as source-sink terms So analytical expressions are exact solution for problem with - these source-sink terms and - first kind boundary conditions with given function values Only good if source-sink terms are small and do not dominate the problem Analytical expressions for „exact“ solution (steady state) Pressure p(x, z) (1 z)(g0 z z z 0 ) x ((3 x s ) (x x s )z) Salt mass fraction 2z 2 2 4.5 2 ( x 3 z) ( x z) x h c( x , z ) 1 e 2 1 Values in example tuned to make sources/sinks small: =20, =12, h=.14, =1, 0=1, =1, x=0.1, z=0.02, xs=1, zs=-0.1 In PDE: n=1, g=1, k=10, =1, Dm=1, /c=0.1, (c)=0+ /c c=1+0.1c Analytical Pressure Salt mass fraction Values between 0 and 1.13 p units Values between 0 and 1 c units Plugged into equations for c and p Source-sink distributions Total mass Salt mass Red: max. input Turquoise: 0 input Red: max. input Light blue: 0 input Blue: output Computed (with 4 grid levels) Pressure Salt mass fraction Values between 0 and 1.13 p units Values between 0 and 1 c units Error Pressure Red : computed value too large by 0.004 % Blue: computed value too small by 0.005 % Salt mass fraction Red : computed value too large by 0.007 % Blue: computed value too small by 0.006 % Experimental benchmark • 3D transient experiment in box with simple boundary and initial conditions • Measurement of concentration distribution in 3D with Nuclear Magnetic Resonance Imaging • Measurement of breakthrough curves Drawback: Test of both model equations and mathematics Way out: Construction of a grid convergent solution inspired by the physical experiments Experimental setup • • • • • • Cube filled with glass beads of diameter 0.7 mm Size of model 20*20*20 cm3 Injection of dense fluid on bottom center hole Application of base flow via top corner holes In unstable case: Injection from below and rotation All parameters measured except transverse dispersivity, diffusion coefficient Experimental setup continued 20 cm 20 cm unstable situation stable situation Salt water Stable situation: Low concentration contrast NMR images of diagonal section: stable situation at low concentration contrast Injection Flushing Equilibration End NMR images of diagonal section: stabel situation at high concentration contrast Injection Equilibration „Entraining“ End Two modes No density contrast Large density contrast Experimental breakthrough curves Low contrast High contrast Comparison of computed and measured breakthrough curves Low contrast High contrast Choice of parameters within intervals given through measurements of those Comparison of concentrations along diagonal section Low contrast case, end of experiment Comparison of concentrations along diagonal section High contrast case, end of experiment Grid convergence: Low contrast case c1i , j 1 c1i , j , c1i 1, j c1i , j 2 2 unit 10-2 Level # grid points 0 8 1 27 2 125 3 729 4 4,913 5 35,937 6 274,625 7 2,146,689 x at level 7: 1.56 mm Grid convergence: Low contrast case Grid convergence: High contrast case c2i , j 1 c2i , j , c2i 1, j c2i , j 2 unit 10-2 Level 0 1 2 3 4 5 6 7 8 # grid points 8 27 125 729 4,913 35,937 274,625 2,146,689 16,974,593 x at level 8: 0.78 mm Grid convergence: High contrast case Error of grid-convergent solutions Low contrast High contrast Unstable case: NMRI data NMRI of vertical and horizontal section: fingering experiment Vertical section Horizontal section Fingering: Experiment vs. Computation Finger velocity: comparison of experiment and computation Causes for poor perfomance • Numerical dispersion smoothes out fingers and eliminates driving force • Initial perturbance not known well enough • Start of fingers on microlevel, not represented by continuum equations Influence of heterogeneity on density flow • Homogeneous Henry problem • Heterogeneous Henry problem Definition of Henry problem: Homogeneous aquifer Hydraulic conductivity 1E-2 m/s, effective diffusion coefficient 1.886E-05 m2/s Boundary conditions: Left fresh water flux given at 6.6E-05 m/s Right hydrostatic salt water, salt mass fraction 0.0357 kg/kg Solution of Henry problem: Homogeneous aquifer Relative concentration contours between 0 (left) and 1 (right) in steps of 0.1 Heterogeneous Henry problem: Permeability distribution Lognormal distribution, exponential autocorrelation Arithmetic mean 1.68E-9 m2, Geometric mean 1.02E-9 m2 Variance of log(k) = 1 Corr. lengths: horizontal 0.05 m, vertical 0.05 m Red: 3.5E-08 m2, Blue: 2.6E-11 m2 Heterogeneous Henry Problem: Concentration distribution Relative concentration contours between 0 (left) and 1 (right) in steps of 0.1 Eff. diffusion coefficient as in homogeneous case Question: Are there equivalent effective parameters to mimick main effect of heterogeneity in a homogeneous model? Heterogeneous Henry Problem Comparison (zero local dispersion) Heterogeneous case Only diffusion Homogeneous case: Permeability equal arithmetic mean of heterogeneous case Only diffusion Homogeneous case: Permeability equal geometric mean of heterogeneous case Only diffusion Heterogeneous Henry Problem Comparison (zero local dispersion) Heterogeneous case only diffusion Homogeneous case: Permeability equal geometric mean of heterogeneous case only diffusion, zero dispersion Homogeneous case: Permeability equal geometric mean of heterogeneous case, diffusion plus macrodispersion after Gelhar&Axness Heterogeneous Henry Problem Comparison (with local dispersion) Heterogeneous case Diffusion + local dispersion Homogeneous case: Permeability geometric mean Diffusion + local dispersion Homogeneous case: Permeability geometric mean Diffusion + macrodispersion after Gelhar&Axness Effective dispersion parameters - Stable situation with flow against direction of density gradient effective longitudinal dispersivity given by Gelhar & Welty, A11(with density gradient) < A11(without density gradient) - Unstable situation with flow in direction of density gradient at sufficiently large Wooding/Raleigh number, dispersivities grow to infinity due to fingers forming - Horizontal flow towards a fixed concentration leads to „boundary layer“. Dispersion at upper right boundary and at stagnation point is „upstream diffusion“: c/c0 = 1 – exp (-x/aL) Upstream diffusion Conclusions • Density flow of increasing importance in groundwater field • Tests for reliability of codes are available • Density flow especially with high contrast numerically demanding: Grid convergence may require millions of nodes • Numerical simulation of fingering instabilities still inadequate • Heterogeneities can be handled by effective media approach in situation without fingering • New numerical methods are in the pipeline