Report

An introduction to PFLOTRAN and its application to CO2 geological storage problems Edinburgh, 10 January 2014 Paolo Orsini PFLOTRAN overview Parallel n-phases n-components reactive flow code for modeling subsurface processes, developed by the cooperation of four US national labs (LANL, PNNL, ORNL, LBNL): Open Source GNU Lesser General Public License (LGPL) Object oriented programming (F95, F2003-2008): easy to extend and to incorporate additional processes Parallel computation based on the PETSC library (ANL lab) Parallel implementation tested on computer architectures with >100k processor cores PFLOTRAN modeling capabilities Solution of mass balance and energy equations that can be coupled sequentially to reactive-transport and quasi-static geo-mechanical models: Single phase variable saturated flow (Richards equation) TH (Thermal Hydrologic) Single phase variable saturated flow with variable density (function of p and T) Immiscible CO2 – brine: non-isothermal two-phase flow Supercritical CO2 – brine: non-isothermal two-phase two-components flow (Variable switch) Development of a black oil model (FVFs) is at planning stage Example of parallel performance on a super computer: Richards Equation (Hanford 300 Area) Variable saturated flow problems 270 M DOF Time[s]: wall clock time per time step Example of parallel performance on a super computer: CO2-Brine CO2-Brine: 25 M cells Yellowstone machine: 8000 cores Flow: 3 DOF Transport: 10 DOF Time[s] for 10 flow step + 14 transport steps PETSc (Portable extensible toolkit for scientific computing) parallel framework - overview Data structure for a parallel (CSR matrices, Blocked CSR matrices, distributed arrays, etc) Non linear solvers (Newton-based methods) Pre-conditioners ( Additive Swartz, ILU, etc..) Time stepper algorithms (Euler, Backward Euler, etc) Krylov Subspace methods (GMRES, CG, CGS, etc) PFLOTRAN Flow diagram Discretisation Space: Control volume (structured and unstructured grids), two point flux formula, MFD under development Time: Flow solvers: implicit Reactive transport Fully implicit Operator splitting (require less memory but also to satisfy the CFL condition for stability Coupling Sequential between flow and reactive transport Domain decomposition Parallelisation based on overlapping-domain decomposition (each processor is assigned to a sub-domain): accumulation terms are easy to compute because are local operations, computation of fluxes require message passing Domain decomposition To evaluate a local function f(x), each process requires the local portion of x and its ghosted part (overlapping part) General n-phase n-components mass balance and energy equations Mass balance equation s X i q X i s D X i Qi t Energy equation sU 1 r crT q H T Qe t q kk P gz W sα phase saturation, η molar density, Xiα molar fraction of component i in phase α, Dα phase diffusivity coefficient, φ porosity, H enthalpy, U internal energy, ρr rock density, κ thermal diffusivity coefficient, kα relative permeability, k saturated media permeability General n-phase n-components mass balance and energy equations: degrees of freedom Gibbs phase rule F Nc 2 N p Open System 2 N p N p Nc Unknowns p, T , s , X i Constraints s 1 X 1 i i i 1 N p N p 1 N c i Degrees of freedom (Ndof) Nc 1 MPHASE - CO2-Brine module: governing equations Two phase [gas, liquid], two components [CO2, H2O]: 3 DOF Two component mass balance: sl l X wl sg g X wg FW QW t sl l X Cl sg g X Cg FC QC t FW ql l X wl q g g X wg l sl Dl X wl g sg Dg X wg FC ql l X Cl q g g X Cg l sl Dl X Cl g sg Dg X Cg MPHASE – CO2-Brine module: auxiliary variables The pure CO2 properties, which depend on P and T, are computed with the Spang & Wagner (1994) EOS. They are tabulated before the computation, and a look-up table is used during the simulation Brine/CO2 mixture density Duan (2008): l P, T , x Solubility of CO2 in brine, Duan and Sun (2003): l g CO CO 2 Solution procedure by variable switch approach l pg , T , sg liquid 2 ph : pl , T , X CO 2 gas 2 ph : g pg , T , X CO pg , T , sg 2 l 2 ph liquid : pg , T , sg pl , T , X CO 2 2 ph gas : g pg , T , sg pg , T , X CO 2 l g CO CO 2 2 l g CO CO 2 sg 0 sg 1 2 2 MPHASE CO2-Brine Module: pressure and temperature range limits The code standard release is limited to Supercritical CO2, however the real limit is the number of phases (CO2 liquid and gas cannot coexist) Geomechanics Model Governing equations (quasi-static equilibrium) b 0 tr 2 P P0 I T T0 I T 1 u x u x 2 u x u p x D x n x t p x N λ and μ are the Lame parameters, related to the Young’s modulus and Poisson ratio. α= coefficient of thermal expansion β=Biot’s coefficient Geomechanics model – Discretisation Equation is solved with the Galerkin finite element method. One-way coupling with the flow solver via pressure and temperature, which are available at the control volume cell centres. The geo-mechanics does not need to be solved at every flow time step. The cell centres are the nodes of the finite element mesh, so there is no need of interpolating P and T. (Voronoi mesh) CV mesh for flow solution FEM mesh geomechanics PFLOTRAN – Pre-Post processors There is no specific pre-processor Geological model and grid generation with external software Several mesh formats: (i) structured with variable spacing [internal generator], (ii) implicit unstructured [list of nodes and connectivity table for hex, wedges, pyramids and tets], (iii) explicit unstructured [list of cell centre volumes and connecting faces, for general polyhedrons] Simulation control parameters, BCs and ICs via flexible text files Post-processor Open source, VisIt, ParaView. Both can post-process remotely, on parallel architectures (auto-reassembling). Commercial software: Tecplot Input deck file to set up the simulation control parametres – organised in cards Grid (internal mesh generator or external mesh) Specify flow mode (the application module: e.g. CO2-brine, Richards) Material properties Capillary & relative permeability functions Regions: interior domain and surfaces Geometry may be specified independent of grid Initial & boundary conditions, source/sinks for flow and transport Coupler: to couple regions with initial and boundary conditions Solvers (direct, Iterative) Time stepping Output Checkpoint/Restart Input deck file example: grid, regions, material, strata Each card containing multiple instructions starts by its key word and ends by “END” or “/” MODE MPHASE GRID TYPE unstructured gridfile.dat ORIGIN 0.d0 0.d0 0.d0 END MATERIAL_PROPERTY soil1 ID 1 POROSITY 0.15d0 TORTUOSITY 1d-1 ROCK_DENSITY 2.65E3 SPECIFIC_HEAT 1E3 THERMAL_CONDUCTIVITY_DRY 2.5 THERMAL_CONDUCTIVITY_WET 2.5 SATURATION_FUNCTION sf2 PERMEABILITY PERM_X 1.d-15 PERM_Y 1.d-15 PERM_Z 1.d-17 / / SATURATION_FUNCTION sf2 PERMEABILITY_FUNCTION_TYPE MUALEM SATURATION_FUNCTION_TYPE VAN_GENUCHTEN RESIDUAL_SATURATION LIQUID_PHASE 0.1 RESIDUAL_SATURATION GAS_PHASE 0.0 LAMBDA 0.762d0 ALPHA 7.5d-4 MAX_CAPILLARY_PRESSURE 1.d6 / REGION reservoir COORDINATES 0. 0. 0. 100. 100. 10. / / / STRATA REGION reservoir MATERIAL soil1 END Input deck file example: flow conditions & initial and boundary condition couplers A line of the input deck file can be commented using “:”, a block using the “skip”/“noskip” key words. The cards can be inserted in any order. MODE MPHASE FLOW_CONDITION initial UNITS Pa,C,M,yr TYPE PRESSURE hydrostatic TEMPERATURE dirichlet CONCENTRATION dirichlet ENTHALPY dirichlet / PRESSURE 2D7 2D7 TEMPERATURE 100 C :TEMPERATURE 75 C CONCENTRATION 1d-6 M ENTHALPY 0.d0 0.d0 END INITIAL_CONDITION FLOW_CONDITION initial REGION all END BOUNDARY_CONDITION west_bound FLOW_CONDITION initial REGION West END skip BOUNDARY_CONDITION west_bound FLOW_CONDITION initial REGION West END Noskip Input deck file example: output, restarts, observations OUTPUT :MASS_BALANCE :TIMES d 0.0 0.1 1.0 TIMES d 0.0 1.0 10. 30. 365.0 730.0 1460.0 FORMAT TECPLOT BLOCK PERMEABILITY POROSITY FORMAT HDF5 PERIODIC_OBSERVATION TIME 50.0 d :PERMEABILITY TIME 1.0 d VELOCITIES / CHECKPOINT 200 RESTART Inj20_pc0-2000.chk REGION well1 COORDINATES 1000.0 1500.0 -1075.5 / END OBSERVATION well1 REGION well1 AT_CELL_CENTER END Geomechanics Model – Demo test case 3D domain with CO2 being injected at the centre of the domain in a deep aquifer formation. Deformation in the domain is considered due to injection. overburden caprock reservoir basement Geomechanics Model – Demo case parameter Problem domain: 2500 x 2500 x 1000 m (x y z) Grid resolution 21 x 21 x 20 for subsurface grid CO2 injection rate: 10 kg/s Young's modulus: 10 Gpa (sandstone) Poisson's ratio: 0.3 Biot Coefficient: 1.0 Coefficient of Thermal Expansion: 10-5 Pa/K Total injection time: 20 y Simulation time: 100 y Displacement in normal directions are set to zero. Top boundary face is free to move Geomechanics Model – Demo case: CO2 saturation Geo-mechanics Model – Demo case: relative vertical displacement Case study: Sleipner – commercial CO2 storage site Reservoir: Utsira formation at a depth of 800-1100 m, average porosity 36%, permeability range from 1000 to 5000 md. Residual gas saturation = 0.21. Many horizontal intra-formation shale layers (0.5 – 2 m thick) that affects the CO2 flow through the reservoir Caprock, shale units with a low permeability of ~ 0.001 md CO2 just above critical conditions on the uppermost layer (Pressure ~80 bars, Temperature ~ 29-33 C) CO2 injected over a 38 m interval of a deviated well at 1012 m depth Injection of about 1Mton of CO2 per year since 1996 Sleipner - Seismic cross section - 2008 Layer 9 Sleipner L9 layer – benchmark released by STATOIL Uppermost point L9 model = -800 m b.s.l. Sea bed ~80 m b.s.l. (T=7 °C). Injection location, (spill/leakage from underneath layer): x~1600m, y~2100m. Injection location MESH and topography (vertical direction out of scale with horizontal direction) Grid: dx=dy=50m, dz~1m (17 layers). Unstructured grid (~310k) prisms. A mesh created by STATOIL for ECLIPSE was converted to the PFLOTRAN format. Parameters controlling the CO2 plume location and distribution Caprock topography Mass rate from the underneath layer (volume rates estimated by structural analysis and seismic measurements Chadwick & Noy 2010) Temperature. Variations close to the CO2 critical temperature value cause large changes in density and viscosity (->mobility) Permeability changes with phase saturation. The relative permeability parameters used in SPE-134891 are adopted (Corey 1958). Capillary pressure has been neglected as suggested in SPE-134891. CO2 properties Computed with the Spang and Wagner EOS implemented in PFLOTRAN: Density limits 500-700 kg/m3 as suggested by Alnes et al (2011) Viscosity doubles with 3 °C temperature reduction Numerical model set up Initial conditions: Hydrostatic pressure: ~ P[8.24, 7.98] MPa Temperature: (a) T=35 °C, rho~500 m3/kg; (b) T= 29 °C, rho~700 m3/kg (Alnes et al 2011) Boundary conditions: (i) top and bottom layer considered perfectly impermeable (replaced by zero flux condition), (ii) side boundaries hydrostatic pressure, (iii) Injection temperature (a) -> T=35 °C, (b) -> T=29 °C Material properties. Rock thermal properties: rho=2600 kg/m3, cp=920 J/kg/°C, ĸ=2.51 W/m/°C (saturated medium). Permeability and porosity assigned cell by cell (Kh~10-12 m2, η~0.35), Kv/Kh=0.1. PFLOTRAN – history matching Top down view of the layer-9 CO2 saturation contour 2002 2004 2006 2008 Monitoring data acquired via repeated seismic and gravimetric surveys by STATOIL ECLIPSE 100 (Oil – gas system) - history matching To conclude PFLOTRAN is distributed by BitBucket, a distributed version control system (DVCS): https://bitbucket.org/pflotran/pflotrandev/wiki/Home PFLOTRAN website: www.pflotran.org User-group forum: google group We are planning to the development of a black oil model within PFLOTRAN: we welcome any suggestion regarding the best type of formulation to implement, and also ideas on how to fund the development Thank you Energy Equation – Thermal conduction multiple continuum model Can be used in the THC and CO2brine flow modules Dual continuum. Primary: fracture (f), ε=Vn/V; secondary matrix (m) Different shape are available for the secondary media (nested sphere, slabs, etc) The solution for the temperature in the secondary media is local (1D) and easy to parallelise sU 1 r crT f t r crTm t Tm f 0 Tm q H f T f Afm fm n PETSC framework – flow of control PFLOTRAN Geo-Mechanical model Geomechanics model – Demo case: Cross section plane; relative displacement vectors at 20 years x rel. displacement z rel. displacement Continuum Damage Formulation (not implemented yet): Empirical model for the propagation of micro-fractures: A.P.S. Selvadurai. Introduces a continuum hydromechanical damage parameter, D. Hydromechanical damage modifies the bulk modulus and permeability. Evolution of damage: Equivalent shear-strain: Where: 2 d dD 1 d d 1 2 d D 1 Dc d eij eij 1 eij ij kk ij 2 1 2 1 ui u j ij 2 x j xi Damage grows where the material deformations are dilatational: tr ij 0 Two primary effects of damage: Increase in hydraulic conductivity: Reduction in bulk modulus: kd 1 d k0 2 d 1 D0 Model is valid up to a critical damage; Dc, beyond which the porous skeleton breaks down and macro-fractures dominate. Representative parameters for sandstone (Selvadurai): 1 2 130 3105 D0 0 DC 0.75 Example: Injection into sandstone at over-pressurisation of 200 m H20 Damage Increase in hydraulic conductivity Large increase in hydraulic conductivity as the injection location is approached. Pre-damaging the injected medium in this way allows a lower injection pressure for a given mass flux. PFLOTRAN Black oil model Development plan Black oil module – 3 phases three components - isothermal 3 Components 3 Phases Oil -> Oil Gas -> Oil Gas -> Gas Water -> Water (S w w ) div( w qw ) Qw t (So o ) div ( o qo ) Qo t (S g g So dg ) div ( g q g dg qo ) Qg Qdg t 200 180 160 140 120 100 80 60 40 20 0 1.4 1.4 1.35 1.3 1.2 1.3 1.1 1.25 1 1.2 0.9 1.15 0.8 0.7 1.1 1.05 0 100 200 300 Bo Vo Vdg RC Vo STC 0.6 muo 0.5 1 400 0 100 P Vdg RS V o STC Bo 200 300 0.4 400 P (bars) ; Vg RC Bg Vg STC RC: reservoir condition, STC: stock tank conditions, Vdg: gas component in the oil phase ; VW RC Bw VW STC muo (cP) Bo (Sm3/m3) Rs(m3/m3) Black oil module – formation volume factor approach Black oil model - development plan Black oil model is valid for dry gas, with a small percentage of volatile component dissolved in the oil phase. The same mathematical formulation can be used for wet and gas condensate with a small oil vapour component Implementation of look-up tables to load the properties required (FVFs, Rs , viscosity) that depend on pressure. The module will be fully integrated into the PFLOTRAN parallel framework and released under the same open source licence Any feedback on what you are not happy with the software you are using at the moment, or any other suggestion is very welcome PFLOTRAN CO2/BRINE Module 2 phases – two components Span and Wagner accuracy Reactive Transport Model Conservation equation for primary species S j s j q s D j Q j jm I m t t m Equilibrium mass action laws for aqueous species Nsec j l C ji Ci l j i 1 Mineral kinetics (transition state theory) n n3,m I m am k0ma exp t f Eama / R aH1,m k0mn exp t f Eamn / R k0mb exp t f Eamb / R aOH tf 1 T 1 298.15 m V m Im t 1 m m k c k0 0 c a am m a0 0 n Black oil module Example of parallel performance on a smaller cluster Transport Equation: 4k cells, 12 species (~50k DOF) Memory contention issues in the first 8 cores, good scalability after Time for each transport step