Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0077 Acta Polytechnica 61(SI):77–88, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague COMPUTATIONAL METHODOLOGY TO ANALYZE THE EFFECT OF MASS TRANSFER RATE ON ATTENUATION OF LEAKED CARBON DIOXIDE IN SHALLOW AQUIFERS Radek Fučíka, ∗, Jakub Solovskýa, Michelle R. Plampinb, Hao Wuc, Jiří Mikyškaa, Tissa H. Illangasekared a Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Mathematics, Trojanova 13, 12000, Praha, Czech Republic b U.S. Geological Survey, Eastern Energy Resources Science Center, 12201 Sunrise Valley Drive, Reston, VA 20192, USA c Virginia Polytechnic Institute and State University, Department of Geosciences, 926 West Campus Drive, Blacksburg, VA 24061, USA d Colorado School of Mines, Center for Experimental Study of Subsurface Environmental Processes, 1500 Illinois St., Golden, CO 80401, USA ∗ corresponding author: fucik@fjfi.cvut.cz Abstract. Exsolution and re-dissolution of CO2 gas within heterogeneous porous media are investigated using experimental data and mathematical modeling. In a set of bench-scale experiments, water saturated with CO2 under a given pressure is injected into a 2-D water-saturated porous media system, causing CO2 gas to exsolve and migrate upwards. A layer of fine sand mimicking a heterogeneity within a shallow aquifer is present in the tank to study accumulation and trapping of exsolved CO2. Then, clean water is injected into the system and the accumulated CO2 dissolves back into the flowing water. Simulated exsolution and dissolution mass transfer processes are studied using both near- equilibrium and kinetic approaches and compared to experimental data under conditions that do and do not include lateral background water flow. The mathematical model is based on the mixed hybrid finite element method that allows for accurate simulation of both advection- and diffusion-dominated processes. Keywords: Compositional flow, two–phase flow, kinetic mass transfer, gas exsolution, gas dissolution. 1. Introduction Geologic carbon sequestration has the potential to significantly reduce greenhouse gas emissions [1], but also poses risks to groundwater resources including mobilization of contaminants in shallow aquifers due to leakage of CO2 from deep storage formations [2]. The extent and severity of the risks depend on complex multiphase flow and transport phenomena that govern the migration of CO2 through the shallow subsurface. A persistent issue with predicting these processes is the general difficulty of understanding CO2 phase change (a.k.a. inter-phase mass transfer) within macroscopic porous media systems, which is important in the case of CO2 due to its high solubility and potential mobility in the gas phase. If a leakage pathway is encountered, stored CO2 that is originally supercritical in a deep geologic stor- age formation may migrate upward due to buoyancy, dissolve into water, exsolve to form a separate gas phase, and eventually re-dissolve into clean water. These interrelated processes are collectively referred to as multiphase evolution, and have recently been stud- ied in various continuum-scale systems. Experimental investigations have identified the factors that control multiphase evolution during one-dimensional (1D) ver- tical flow [3–5] and have qualitatively addressed the various transport phenomena during two-dimensional (2D) flow [6]. Numerical modeling has been used to show the effects of 1D flow rate [7] and quasi-2D flow paths [8] on multiphase CO2 evolution processes. How- ever, numerical models have not yet been able to fully explain all of the observations from the experimental studies, particularly of those that occur during 2D flow under the influence of background water flow. The purpose of this present study is to develop a numerical modeling methodology to assess the factors that control inter-phase CO2 mass transfer during migration through heterogeneous 2D porous media systems in the shallow subsurface where CO2 can only exist in the gaseous and dissolved phases. Specifically, we seek to test whether the local equilibrium assump- tion for mass transfer adequately explains multiphase evolution, and to investigate the effect of lateral water flow on heterogeneity-enhanced gas phase accumula- tion. Our approach is to expand the numerical model developed in [9] to include kinetic mass transfer, and then to compare the results of the simulations per- formed by the model to the laboratory data from experiments that build upon methods developed in 77 https://doi.org/10.14311/AP.2021.61.0077 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica [10]. This methodology forms the basis for another article [11] that investigates multiphase CO2 plume dynamics in relatively large-scale synthetic aquifer sys- tems. These investigations improve our understanding of the underlying processes that together lead to at- tenuation of CO2 within groundwater systems [6]. The main objective of this work is to study the exsolution and dissolution of CO2 on a small scale using both experiments and numerical modeling. Mo- tivated by fundamental differences in mass transfer in various scenarios studied in [11], our aim is to de- termine whether a kinetic or a simplified equilibrium mass transfer model is needed to reproduce the results observed in the experiments and how the results are affected by changes in the flow field enforced in the experiments. 2. Experiment The experimental methodology is detailed by [10], so it will be only briefly described here. To assess the effects of background (lateral) water flow on multiphase CO2 migration, two different bench-scale experiments were performed. In each experiment, CO2-saturated water was injected into the bottom of a tank filled with water-saturated porous media until gas phase formed and accumulated to some steady value. Then, clean water was injected until all of the CO2 dissolved and migrated out of the tank. The test system was initially reported by [10], but this study expands upon that methodology by in- corporating different experimental conditions. Most importantly, one experiment was performed with wa- ter flowing from right to left across the tank. As is shown in Figure 1, a block of low-permeability sand was incorporated into the middle of the tank above the injection ports, and saturation sensors were installed below this block of fine sand where gas phase CO2 was expected to accumulate after exsolving from the injected CO2-saturated water. Material properties of the three porous media used are given in Table 3. The saturation pressure, as defined by [3], was 10 kPa for the first experiment and 15 kPa for the second, which led to dissolved CO2 concentrations that were sufficiently high to cause exsolution in the porous media. In the first experiment, the constant head devices connected to the regions of Granusil #8 (which were installed to distribute the head evenly across the verti- cal boundaries of the main sand pack) were positioned at equal elevations. This led to negligible background water flow across the tank, and this case is therefore re- ferred to as the static case. In the second experiment, the left hand constant head device was positioned at an elevation below that of the right-hand one, thus establishing background flow. In the second case, wa- ter inflow was supplied to the right-hand constant head device via a peristaltic pump, while the outflow from the left-hand constant head device was routed into a container placed on a computer-interfaced elec- tronic scale. The scale and the saturation sensors were configured to automatically take measurements at 1-minute intervals. The resulting data from the experiments were com- pared against a numerical model of two-phase flow in porous media, which extends upon that of [9]. A major addition to the model is its capability to ac- count for non-equilibrium CO2 mass transfer between aqueous and gas phases. 3. Mathematical model In this section, the mathematical model that describes the two-phase compositional flow in porous media and incorporates the phenomena studied in this work is summarized. 3.1. Two-phase flow in porous media The governing equations for the two-phase flow in porous media are based on [12–14]. The quantities corresponding to the liquid (wetting) and gas (non- wetting) phases are denoted by indices ` and g, re- spectively. The mass balance equations for the incompressible liquid and gas phases are given by ρ`φ ∂S` ∂t + ρ`∇·~v` = f`, (1a) and ρgφ ∂Sg ∂t + ρg∇·~vg = fg, (1b) respectively, where φ [−] is the material porosity and Sα [−], ρα [kg m−3], ~vα [m s−1], fα [kg m−3 s−1] are the α-phase saturation, density, velocity, and the sink or source term. The velocity ~vα is given by Darcy’s law ~vα = −λαK(∇pα −ρα~g), (2) where ~g [m s−2] is the gravitational acceleration vec- tor, K [m2] is the intrinsic permeability, pα [Pa] is the α-phase pressure, λα = krα/µα [Pa −1 s−1] denotes the mobility of phase α, where µα [Pa s] is the dy- namic viscosity, and krα(Sα) [−] denotes the relative permeability. The difference between the wetting and non-wetting phase pressures is defined as the capillary pressure pc = pg −p` and the Brooks and Corey model [15] is used in the form pc(S`) = pd(Se` ) − 1 λ , (3) where pd [Pa] is the entry pressure, λ [−] is the pore size distribution index, and Seα [−] denotes the effec- tive saturation defined by Seα = Sα −Sr,α 1 −Sr,g −Sr,` , (4) where Sr,α [−] is the residual saturation of phase α, α ∈{`,g}. 78 vol. 61 Special Issue/2021 Computational Methodology of Mass Transfer on Attenuation of CO2 Water outflow (to scale) Water outflow (to scale) CO2 gas outflow (flow meter) 5 .8 c m Granusil #110 & #250 (2:1) Granusil #20/#30 Granusil #8 (gravel)Saturation sensor CO2-saturated water injection port Temperature sensor 3.81 cm 3.81 cm 25.4 cm 3.81 cm 3.81 cm 1 5 .2 4 c m 1 5 .2 4 c m 1 cm Headspace 1 cm 2.5 cm 10.2 cm 2.5 cm 1.5 cm port C port B port A port D port E port F Figure 1. Setup for the small tank experiments, adapted from [10]. The Brooks and Corey model parameters are also used in the Burdine’s model for the relative perme- ability functions kr,` and kr,g, [16], in the form kr,`(S`) = (Se` ) 2+3λ λ , (5) kr,g(Sg) = (Seg) 2 ( 1 − (1 −Seg) 2+λ λ ) . (6) Based on [7, 11, 17], however, instead of kr,g, a mod- ified formula for the gas phase relative permeability function is used in Eq. (2) in the form k̃r,g(Sg) = { 0, if Sg < Sc, kr,g( Sg−Sc 1−Sc ), otherwise, (7) where Sc [−] denotes the critical gas saturation. 3.2. Component transport The gas phase is considered as a single component (pure CO2), whereas the liquid phase is assumed to be a two component mixture (water and CO2). The compositional balance equation for CO2 dissolved in the liquid phase is included into (1) and (2) as φρ` ∂(S`X) ∂t + ρ`∇· (X~v` − τφS`D`∇X) = fX, (8) [14], where X [−] is the mass fraction of CO2, fX [kg m−3 s−1] is the sink or source term, D` [m2 s−1] is the free molecular diffusion of CO2 in water, D` = 1.92 · 10−9 m2 s−1, and τ` [−] is the tortuosity given by τ` = φ1/3S 7/3 ` based on [18]. 3.3. Kinetic mass transfer Based on [19], the kinetic mass transfer of CO2 be- tween both phases (i.e., the dissolution and exsolution processes) is represented by −fg = f` = fX = k(Cs −Xρ`), (9) where k [s−1] is the lumped mass transfer rate coeffi- cient and Cs [kg m−3] is the saturated CO2 concentra- tion in water at the relevant pressure and temperature given by Henry’s law in the form Cs = pg KH Mg, (10) 79 R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica where Mg [ kg mol−1 ] is the molar mass of CO2, Mg = 44.01 g mol−1, and KH [ Pa mol−1 m3 ] is Henry’s constant for which the Van’t Hoff equation is employed in the form KH = KH,refe −C ( 1 T − 1 Tref ) (11) where T [K] is the temperature, KH,ref is the value of Henry’s constant at a reference temperature Tref [K], and C [K] is the gas-specific constant, i.e., KH,ref = 2979.97 [ Pa mol−1 m3 ] , Tref = 298.15 K, and C = 2400 K [20]. The average temperatures in the static case and background flow experiments were 37◦C and 26◦C, respectively. 4. Numerical method and implementation The governing equations are solved using a general numerical solver based on the mixed-hybrid finite element method described in [9]. The mixed-hybrid fi- nite element method combines velocity discretizations in the lowest order Raviart–Thomas–Nédélec space with piecewise constant approximations for the scalar variables. The numerical method can be used for ac- curate simulation of degenerate diffusion or advection- dominated problems like the one discussed here. 4.1. Numerical method The numerical method solves a system of n partial differential equations in the coefficient form in a d- dimensional polygonal domain Ω ⊂ Rd and a time interval [0, tfin]: n∑ j=1 Ni,j ∂Zj ∂t + n∑ j=1 ~ui,j ·∇Zj+ ∇· (mi(~qi + ~wi)) = fi, (12a) with ~qi = − n∑ j=1 Di,j∇Zj, (12b) where Zj = Zj(t,~x), j = 1, 2, . . . ,n, represent the unknown variables, ~x ∈ Ω, t ∈ [0, tfin]. Eq. (12) is fur- ther supplemented with Dirichlet or Neumann bound- ary conditions, or its combination on different parts of the boundary [9]. Here we restrict ourselves only the case d = 2 which is relevant for the problems presented in the paper. In brevity, the mixed hybrid finite element discretization approximates the unknown scalar functions Zj in the space of piecewise constant functions, i.e., Zj ≈ ∑ A∈Ah Zj,AϕA (13) and the vector function ~qi from Eq. (12b) in the lowest order Raviart–Thomas–Nédélec space RTN0 as ~qi ≈ ∑ A∈Ah χA ∑ E∈EA qi,A,E~ωA,E, (14) where Ah denotes the set of triangles discretizing the computational domain Ω, h is the mesh size defined as the largest ball diameter circumscribed around elements in Ah, EA is the set of all sides of an element A ∈Ah, ϕA are the piecewise constant basis functions defined by ϕA = χA, where χA is the characteristic function of A defined by χA(~x) = { 1 ~x ∈ A, 0 ~x /∈ A, (15) and ~ωA,E are the vector basis functions of RTN0(A), [9]. Eq. (12b) and the discretization given by Eq. (14) allow to express the coefficients qi,A,E as qi,A,E = ∑ j∈σi,A ( bi,j,A,EZj,A − ∑ F∈EA bi,A,E,FZj,F ) , (16) where Zj,E are traces of Zj on side E ∈ EA, σi,A ⊆ {1, 2, . . . ,n} denotes the set of all indices j for which Di,j is non-zero on element A and the definition of coefficients bi,j,A,E and bi,j,A,E,F are given in [9]. The approximation of ~qi in RTN0 given by Eqs. (14) and (16) together with the piecewise constant approx- imation of Eq. (12a) allow to express the cell-averages Zj,A as a local linear combination of Zj,E, E ∈ EA, on each element A ∈ Ah, see [9] for details. Hence, qi,A.E can be expressed solely as linear function of traces Zj,E, ∀E ∈EA. Across all interior sides E of the triangulation Ah, the balance of conservative fluxes is given by∑ {A:E∈EA} mi,E(~qi,A,E + ~wi,A,E) = 0, (17) where a single (upwinded) value mi,E is assumed at E, see [9], and, thus, it can be eliminated from Eq. (17) if it has a non-zero value to produce∑ {A:E∈EA} ~qi,A,E + ~wi,A,E = 0. (18) When mi,E = 0 in Eq. (17) for some side E, Eq. (18) is also considered to assure that the resulting system of linear equations, given in the vector form as Mk ~Z`+1 = ~b`, (19) is non-singular. In Eq. (19), M is a sparse, positive definite matrix [9] and ~b is the right-hand-side, both evaluated on the previous time level `, and ~Z is the vector containing side-traces Zj,E on the next time level ` + 1 for all j = 1, 2, . . . ,n and all sides E in Ah. 4.2. Problem formulation The system of governing equations given by Eqs. (1), (2), (8), and (9) is represented by (12) using n = 3, 80 vol. 61 Special Issue/2021 Computational Methodology of Mass Transfer on Attenuation of CO2 Parameter Units liquid gas H2O CO2 Density ρ [kg m−3] 997.78 1.98 Dyn. viscosity µ [10−5Pa s−1] 87.2 1.48 Table 1. Properties of fluids. Units Mesh Mesh Mesh ∆1 ∆2 ∆3 Elements 5 494 18 587 72 949 Mesh size h [mm] 5.55 2.78 1.43 Table 2. Properties of meshes. d = 2, Z1 = pc, Z2 = pg, Z3 = X, and ( Ni,j ) i,j∈3̂ =  −φρ` dS`dpc 0 0−φρg dS`dpc 0 0 0 0 φS`ρ`   , (20) ( ~ui,j ) i,j∈3̂ =  ~0 ~0 ~0~0 ~0 ~0 ~0 ~0 ρ`~v`   , (21) ( mi ) i∈3̂ =  ρ` λ` λ`+λg ρg λg λ`+λg ρ`   , (22) ( Di,j ) i,j∈3̂ =  (λ` + λg)K −λ` + λgK 00 (λ` + λg)K 0 0 0 τφS`D`   , (23) ( ~wi ) i∈3̂ =  −(λ` + λg)ρ`K~gλ` + λgρgK~g ~0   , (24) ( fi ) i∈3̂ =   −f`fg f` −Xf`   , (25) where 3̂ = {1, 2, 3}. The computational domain Ω depicted in Figure 2 is discretized using three gradually refined triangular meshes generated by gmsh [21]. The mesh properties are listed in Table 2 and in Figure 2, the coarsest mesh is shown. The numerical method is implemented in an in- house computer code using C++. The applicability of the numerical method for heterogeneous porous media is further discussed in [22] together with parallel implementations of the method on GPU [9] or on CPU using MPI [23]. 5. Results and discussion 5.1. Computational study After the physical domain and boundaries were set up to mimic the experiments, several simulations were performed for each experimental case (i.e., the static Γ1 Γ2 ΓW ΓE ΓW ΓE ΓW ΓE ΓS ΓN Γ3 Γ4 Γ6 Γ5 5 cm 0.8 cm 9.2 cm 0.8 cm 9.92 cm0.4 cm20.12 cm 0.4 cm Figure 2. Computational domain with coarse mesh based on Figure 1. one and the one with background flow), using sev- eral different values for k. The boundary conditions on the model included hydrostatic pressure in the inflow/outflow ports on the left and right-hand side of the tank (see Figure 1). For the experiment with back- ground flow, the pressure on the right-hand side was increased to enforce a lateral pressure gradient and thus background flow through the tank. The pressure difference across the tank was fitted to match the mea- sured outflow rates from the experiment. Injection of CO2-saturated and clean water was represented by a Neumann boundary condition at the injection port. No flow boundary conditions were prescribed along the remaining parts of the computational do- main boundary. Several additional simulations were performed to in- vestigate the effect of different background water flow rates on gas accumulation below the heterogeneity. All of these simulations were performed with a high (near-equilibrium) mass transfer rate of k = 1 s−1. The simulation representing the experiment with back- ground flow was conducted with a pressure difference of 5 Pa, while the other simulations were performed with pressure differences of 2, 4, 6, 7, 8, and 9 Pa. 5.2. Numerical convergence The mesh effects on the numerical solution are investi- gated using three gradually refined triangular meshes described in Table 2. The gas saturation profiles at Ports A, B, and C for k = 0.5 s−1 and Sc = 0.2 for both static and background flow cases are compared in Figure 3. The amount of gas simulated present in ports A, B, and C is slightly overestimated for coarser meshes with respect to finer meshes, however, these results are expected due to the first order of convergence of the numerical method as reported in [9]. 5.3. Discussion In this section, the comparison between the experi- mental data and the numerical results is presented 81 R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica Parameter Units Granusil Granusil 2:1 mixture of Granusil #8 #20/30 #110 and #250 Porosity φ [-] 0.4 0.32 0.35 Intrinsic permeability K [m2] 1×10−9 2.30×10−10 6.36×10−14 Pore size distribution index λ [-] 4.275 7.33 5.35 Entry pressure pd [Pa] 600 1200 8027 Residual liquid phase saturation Sr,` [-] 0.084 0.084 0.17 Residual gas phase saturation Sr,g [-] 0 0 0 Table 3. Material properties. 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 60 70 Time [h] f) Port C, background flow 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 Time [h] e) Port C, static case 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 60 70 Time [h] d) Port B, background flow 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 Time [h] c) Port B, static case 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 60 70 Time [h] b) Port A, background flow 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 10 20 30 40 50 Time [h] a) Port A, static case Mesh ∆1 Mesh ∆2 Mesh ∆3 Figure 3. Illustration of mesh resolution effects on the numerical solution of the gas saturation evolution for k = 0.5 s−1 and Sc = 0.2 at the sampling ports A (a b), B (c, d), and C (e, f) for both static (a, c, e) and background flow cases (b, d, f). The mesh properties are given in Table 2. 82 vol. 61 Special Issue/2021 Computational Methodology of Mass Transfer on Attenuation of CO2 0 0.05 0.1 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 50 60 Time [h] f) Port C, background flow 0 0.05 0.1 0.15 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 Time [h] e) Port C, static case 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 50 60 Time [h] d) Port B, background flow 0 0.05 0.1 0.15 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 Time [h] c) Port B, static case 0 0.05 0.1 0.15 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 50 60 Time [h] b) Port A, background flow 0 0.05 0.1 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 Time [h] a) Port A, static case Model, k = 0.01 s−1 Model, k = 0.05 s−1 Model, k = 0.1 s−1 Model, k = 0.5 s−1 Model, k = 1 s−1 Model, k = 2 s−1 Experiment Figure 4. Results of the experiments and numerical simulations for (a) the static case and (b) the case with background flow. 83 R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica 0 0.1 0.2 G a s sa tu ra ti o n S g [− ] 0 10 20 30 40 50 60 Time [h] Model, 2 Pa Model, 4 Pa Model, 5 Pa Model, 6 Pa Model, 7 Pa Model, 8 Pa Port A Figure 5. Results of the simulations performed with various pressure differences applied across the tank. The experiment with background flow corresponds to the simulation with a pressure difference of 5 Pa. and the mass transfer and accumulation processes are discussed for both static and background flow cases. Figure 4 shows the results of gas saturation from the experiments and simulations for the static case and the case with background flow. The experimental data shown in (a) is from port B (directly above the injection port, as CO2-saturated water was injected into the left injection port in this case), while that in (b) is from port A (also above the injection port, as CO2-saturated water was injected into the right injection port in this case). Both experiments show the accumulation of gas phase CO2, then re-dissolution during the injection of clean water after the CO2 saturations had stabilized. The simulation results for various choices of k show that the best match with experimental data is achieved for rather large values of the mass transfer coefficient and that the simulated results are almost the same for those greater than a certain threshold value of 0.1 to 0.5 s−1. This indicates that the mass transfer pro- cesses are at near-equilibrium and that the background flow is not strong enough to substantially affect the CO2 exsolution and dissolution rates in these cases. On the other hand, the numerical results show that the presence of the background flow slows down the accumulation dynamics below the fine layer as indi- cated by different slopes of the saturation curves in Figure 4c compared to Figure 4b. Note that because different injection ports were used in each experiment, the corresponding ports are shifted by one, i.e., re- sults from port B for the static case are compared to the results from port A for the background flow case (see the gas saturation and CO2 mass fraction spatial evolution shown in Figures 6 and 7). Further investigation of the background flow effect on the exsolution process can be done by analyzing the results of the simulations performed with different background flow rates, as shown in Figure 5. The re- sults show a wide range of different gas accumulation behaviors, with higher flow rates leading to lower and slower exsolution due to different gas phase distribu- tion around the injection ports. This indicates that the flow field significantly affects the multiphase CO2 evolution, however, the experimental data for these scenarios would be needed to quantify the impact. While the analysis above indicated that near- equilibrium behavior occurred in the experiments, further experimental data are needed to determine generally in which cases the equilibrium mass transfer simplification is valid. Both experiments considered in this work are on a small scale and the measurement ports are close to the injection. Also, the difference between the static and background flow cases is not very substantial. Addressing these remaining knowl- edge gaps is one of the primary goals in [11], where non-equilibrium mass transfer behavior is observed in a large laboratory-scale system with background flow, and fundamentally different results are observed as compared with predominantly one-dimensional flows. This indicates that, for the decision between the ki- netic and equilibrium models, flow field and scale both need to be taken into account. Additionally, the flow field and scale are not the only factors controlling the multiphase evolution, as also indicated in [11], which showed that temperature fluctuations can also significantly affect mass transfer processes. 6. Conclusions To assess the potential risks of groundwater contam- ination via leakage of stored CO2 from deep geo- logic carbon sequestration sites, we must be able to predict multiphase CO2 transport through shallow aquifers. To aid in this effort, an innovative numeri- cal model was applied and tested against data from well-controlled, small-scale laboratory experiments. Beyond conventional continuum-based two-phase flow 84 vol. 61 Special Issue/2021 Computational Methodology of Mass Transfer on Attenuation of CO2 Gas saturation Sg, static case time 1 h Mass fraction X, static case time 1 h Gas saturation Sg, static case time 10 h Mass fraction X, static case time 10 h Gas saturation Sg, static case time 20 h Mass fraction X, static case time 20 h Figure 6. Spatial distribution of Sg (left) and X (right) for the static case at 1 h, 10 h, and 20 h (from top to bottom), computed using k = 0.5 s−1 and Sc = 0.2. 85 R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica Gas saturation Sg, background flow time 1 h Mass fraction X, background flow time 1 h Gas saturation Sg, background flow time 10 h Mass fraction X, background flow time 10 h Gas saturation Sg, background flow time 20 h Mass fraction X, background flow time 20 h Figure 7. Spatial distribution of Sg (left) and X (right) for the background flow case at 1 h, 10 h, and 20 h (from top to bottom), computed using k = 0.5 s−1 and Sc = 0.2. 86 vol. 61 Special Issue/2021 Computational Methodology of Mass Transfer on Attenuation of CO2 in porous media, the model incorporated the ability to simulate non-equilibrium CO2 mass transfer between the aqueous and gaseous phases during exsolution and dissolution processes. The experimental data used in this work clearly showed equilibrium or near-equilibrium behavior for both steady and background flow cases. These find- ings are in contrast with the non-equilibrium behavior observed on a larger scale for the scenarios with back- ground flow. This indicates that the mass transfer for the similar scenarios is not controlled by the flow field only. Further investigation of the mass transfer is needed to understand for which scenarios the kinetic model is needed and when the equilibrium assump- tion is sufficient. The computational methodology developed here could be used to assess both equilib- rium and non-equilibrium processes in these future scenarios. 7. Acknowledgments The work reported in this paper was supported by the Ministry of Education, Youth and Sports of the Czech Republic under the OP RDE grant number CZ.02.1.01/0.0/0.0/16_019/0000778 Centre for Ad- vanced Applied Sciences and Inter-Excellence grant no. LTAUSA19021. Any use of trade, firm, or product names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government. References [1] S. Pacala, R. Socolow. Stabilization wedges: solving the climate problem for the next 50 years with current technologies. Science 305(5686):968–972, 2004. doi:10.1126/science.1100103. [2] J. A. Apps, L. Zheng, Y. Zhang, et al. Evaluation of potential changes in groundwater quality in response to CO2 leakage from deep geologic storage. Transport in Porous Media 82(1):215–246, 2010. doi:10.1007/s11242-009-9509-8. [3] M. R. Plampin, T. H. Illangasekare, T. Sakaki, R. J. Pawar. Experimental study of gas evolution in heterogeneous shallow subsurface formations during leakage of stored CO2. International Journal of Greenhouse Gas Control 22:47–62, 2014. doi:10.1016/j.ijggc.2013.12.020. [4] M. R. Plampin, R. N. Lassen, T. Sakaki, et al. Heterogeneity-enhanced gas phase formation in shallow aquifers during leakage of co2-saturated water from geologic sequestration sites. Water Resources Research 50(12):9251–9266, 2014. doi:10.1002/2014WR015715. [5] T. Sakaki, M. R. Plampin, R. Pawar, et al. What controls carbon dioxide gas phase evolution in the subsurface? experimental observations in a 4.5 m-long column under different heterogeneity conditions. International Journal of Greenhouse Gas Control 17:66–77, 2013. doi:10.1016/j.ijggc.2013.03.025. [6] M. R. Plampin, M. L. Porter, R. J. Pawar, T. H. Illangasekare. Intermediate-scale experimental study to improve fundamental understanding of attenuation capacity for leaking CO2 in heterogeneous shallow aquifers. Water Resources Research 53(12):10121–10138, 2017. doi:10.1002/2016WR020142. [7] M. L. Porter, M. Plampin, R. Pawar, T. Illangasekare. CO2 leakage in shallow aquifers: A benchmark modeling study of CO2 gas evolution in heterogeneous porous media. International Journal of Greenhouse Gas Control 39:51–61, 2015. doi:10.1016/j.ijggc.2015.04.017. [8] R. N. Lassen, M. R. Plampin, T. Sakaki, et al. Effects of geologic heterogeneity on migration of gaseous CO2 using laboratory and modeling investigations. International Journal of Greenhouse Gas Control 43:213–224, 2015. doi:10.1016/j.ijggc.2015.10.015. [9] R. Fučík, J. Klinkovský, J. Solovský, et al. Multidimensional mixed–hybrid finite element method for compositional two-phase flow in heterogeneous porous media and its parallel implementation on gpu. Computer Physics Communications 238:165–180, 2019. doi:10.1016/j.cpc.2018.12.004. [10] M. R. Plampin, M. Porter, R. Pawar, T. H. Illangasekare. Multi-scale experimentation and numerical modeling for process understanding of co2 attenuation in the shallow subsurface. Energy Procedia 63:4824–4833, 2014. doi:10.1016/j.egypro.2014.11.513. [11] J. Solovský, R. Fučík, M. R. Plampin, et al. Dimensional effects of inter-phase mass transfer on attenuation of structurally trapped gaseous carbon dioxide in shallow aquifers. Journal of Computational Physics 40, 2020. doi:10.1016/j.jcp.2019.109178. [12] P. Bastian. Numerical Computation of Multiphase Flows in Porous Media. Habilitation thesis, Kiel University, 2000. [13] R. Helmig. Multiphase Flow and Transport Processes in the Subsurface, A contribution to the Modelling of Hydrosystems. Springer, 1997. [14] K. Mosthaf, K. Baber, B. Flemisch, et al. A coupling concept for two-phase compositional porous-medium and single-phase compositional free flow. Water Resources Research 47(10), 2011. doi:10.1029/2011WR010685. [15] R. Brooks, A. Corey. Hydraulic Properties of Porous Media. Colorado State University, Hydrology Paper 3:27, 1964. [16] N. Burdine. Relative permeability calculations from pore size distribution data. Journal of Petroleum Technology 5:71–78, 1953. doi:10.2118/225-G. [17] I. Tsimpanogiannis, Y. C. Yortsos. The critical gas saturation in a porous medium in the presence of gravity. J Colloid Interface Sci 270:388–395, 2014. doi:10.1016/j.jcis.2003.09.036. [18] R. J. Millington, J. P. Quirk. Permeability of porous solids. Transactions of the Faraday Society 57:1200–1207, 1961. [19] J. Niessner, S. M. Hassanizadeh. Modeling kinetic interphase mass transfer for two-phase flow in porous media including fluid-fluid interfacial area. Transport in Porous Media 80:329–344, 2009. doi:10.1007/s11242-009-9358-5. 87 http://dx.doi.org/10.1126/science.1100103 http://dx.doi.org/10.1007/s11242-009-9509-8 http://dx.doi.org/10.1016/j.ijggc.2013.12.020 http://dx.doi.org/10.1002/2014WR015715 http://dx.doi.org/10.1016/j.ijggc.2013.03.025 http://dx.doi.org/10.1002/2016WR020142 http://dx.doi.org/10.1016/j.ijggc.2015.04.017 http://dx.doi.org/10.1016/j.ijggc.2015.10.015 http://dx.doi.org/10.1016/j.cpc.2018.12.004 http://dx.doi.org/10.1016/j.egypro.2014.11.513 http://dx.doi.org/10.1016/j.jcp.2019.109178 http://dx.doi.org/10.1029/2011WR010685 http://dx.doi.org/10.2118/225-G http://dx.doi.org/10.1016/j.jcis.2003.09.036 http://dx.doi.org/10.1007/s11242-009-9358-5 R. Fučík, J. Solovský, M. R. Plampin et al. Acta Polytechnica [20] R. Sander. Compilation of Henry’s Law Constants for Inorganic and Organic Species of Potential Importance in Environmental Chemistry. Max-Planck Inst of Chem, Mainz, Germany 1999. [21] C. Geuzaine, J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering 79(11):1309–1331, 2009. doi:10.1002/nme.2579. [22] J. Solovský, R. Fučík. Mass Lumping for MHFEM in Two Phase Flow Problems in Porous Media. In F. A. Radu, K. Kumar, I. Berre, et al. (eds.), Numerical Mathematics and Advanced Applications ENUMATH 2017, pp. 635–643. Springer International Publishing, 2019. doi:10.1007/978-3-319-96415-7_58. [23] J. Solovský, R. Fučík. A Parallel mixed-hybrid finite element method for two phase flow problems in porous media using MPI. Computer Methods in Materials Science 17:84–93, 2017. 88 http://dx.doi.org/10.1002/nme.2579 http://dx.doi.org/10.1007/978-3-319-96415-7_58 Acta Polytechnica 61(SI):77–88, 2021 1 Introduction 2 Experiment 3 Mathematical model 3.1 Two-phase flow in porous media 3.2 Component transport 3.3 Kinetic mass transfer 4 Numerical method and implementation 4.1 Numerical method 4.2 Problem formulation 5 Results and discussion 5.1 Computational study 5.2 Numerical convergence 5.3 Discussion 6 Conclusions 7 Acknowledgments References