An Adaptive Hybrid Vertical Equilibrium/Full‐Dimensional Model for Compositional Multiphase Flow

Efficient compositional models are required to simulate underground gas storage in porous formations where, for example, gas quality (such as purity) and loss of gas due to dissolution are of interest. We first extend the concept of vertical equilibrium (VE) to compositional flow, and derive a compositional VE model by vertical integration. Second, we present a hybrid model that couples the efficient compositional VE model to a compositional full‐dimensional model. Subdomains, where the compositional VE model is valid, are identified during simulation based on a VE criterion that compares the vertical profiles of relative permeability at equilibrium to the ones simulated by the full‐dimensional model. We demonstrate the applicability of the hybrid model by simulating hydrogen storage in a radially symmetric, heterogeneous porous aquifer. The hybrid model shows excellent adaptivity over space and time for different permeability values in the heterogeneous region, and compares well to the full‐dimensional model while being computationally efficient, resulting in a runtime of roughly one‐third of the full‐dimensional model. Based on the results, we assume that for larger simulation scales, the efficiency of this new model will increase even more.

BECKER ET AL. 10.1029/2021WR030990 2 of 19 on gas plume migration and long-term trapping. Compositional effects also have been developed and tested by Nilsen et al. (2016) for a sharp-interface VE model considering residual CO 2 and dissolution of CO 2 in brine only. Two models were developed: one with instant dissolution of CO 2 into brine along the vertical height of the cell, and a second one with a (constant) rate-dependent dissolution of CO 2 into brine (similar to Gasda et al., 2011b). Results with both developed models show that there is a significant total amount of gas dissolving into brine, which in turn has a strong influence on gas plume development.
The model with a constant dissolution rate was extended by Møyner et al. (2020), with a numerical upscaling to account for the transient dissolution of gas components below the gas-phase plume. For this purpose, lookup tables are generated based on pre-computed one-dimensional simulations that link the vertical height of the column, coarse-scale pressure, temperature, time since onset of dissolution, and total mass of components to the coarse-scale deviation from chemical equilibrium. Furthermore, this framework allows for coupling of a full-dimensional model to a compositional VE model, while subdomains need to be assigned a priori, for example, near the injection region. The model developed by Møyner et al. (2020) proved to be efficient as its results compared well against those from full-dimensional simulations, especially regarding the accumulated mass of gas components in brine and gas phase.
Applicability of VE models in general depends on length and time scale (Bandilla et al., 2019;Court et al., 2012), and on locally relevant processes such as, for example, vertical flow due to heterogeneity (Bandilla et al., 2017). As a result, VE may not be valid in the entire domain during the entire simulation time. VE of phases develops with time due to the density difference between the developing gas phase and the resident brine. In some regions such as close to the injection or at heterogeneities, VE of phases may be reached much later or not at all during simulation time. Chemical equilibrium can be considered to be valid almost instantaneously inside the gas-phase plume, whereas chemical equilibrium along the entire depth of the aquifer is reached only after much larger time scales. Therefore, a model that adaptively matches model complexity to system complexity is desirable.
The main contributions of the present work are: (a) The introduction of a hybrid compositional model that adaptively couples a compositional full-dimensional model in subregions with vertical dynamics to a vertically integrated compositional model everywhere else in the domain and (b) the development and analysis of a local criterion to identify where the compositional VE assumption is valid in the domain, including extraction and hysteretic effects due to residual saturation. The advantage of this new model concept is that it combines both accuracy of the full-dimensional model and efficiency of the VE model in one framework. To this end, the compositional VE model is coupled to the compositional full-dimensional model by formulating the volume change due to fluxes across subdomain interfaces. As the hybrid model runs, VE subdomains are identified by the local criterion and the models are assigned adaptively to these regions. We use a test case of hydrogen storage in a heterogeneous, idealized dome-shaped aquifer, and demonstrate efficiency and accuracy of the compositional VE model and the compositional hybrid model. Furthermore, we vary the permeability of the heterogeneous regions and show how the hybrid model adapts to different local segregation times.
This study is structured as follows. In Section 2, the VE assumption is extended to compositional multiphase flow, and model concepts for compositional VE are discussed. In Section 3, a sequential compositional full-dimensional model is presented and integrated in the vertical direction to derive the compositional VE model. The coupling strategy is presented in Section 4, with a formulation of the volume change due to fluxes across subdomain interfaces and an algorithm that solves the coupled system. Following that, a criterion for compositional VE is developed and analyzed, and the adaptive solution algorithm is presented in Section 5. Finally, Section 6 demonstrates the applicability and efficiency of the developed compositional VE model and the compositional hybrid model on a scenario of hydrogen storage in a porous aquifer.

Vertical Equilibrium for a Compositional Two-Phase System
Three different states of VE can be identified for compositional flow in the context of gas storage: 1. Phase equilibrium: gas and brine phase are in equilibrium, resulting in hydrostatic pressure distribution in the vertical direction.
BECKER ET AL.

10.1029/2021WR030990
3 of 19 2. Chemical equilibrium within the gas-phase plume: everywhere where the gas phase is present, there is chemical equilibrium between the gas and brine phase contained in this region. Chemical equilibrium means that, locally, the chemical potential between the present phases is equal. 3. Chemical equilibrium over the entire depth of the aquifer. This includes chemical equilibrium within the gas-phase plume as described above, as well as a fully gas-saturated brine phase below the gas-phase plume.
We expect the time scales for these three equilibrium states to differ substantially, with phase equilibrium being reached first, followed almost instantaneously by chemical equilibrium within the gas-phase plume. The time scale associated with chemical equilibrium over the entire depth of the aquifer would then be expected to be much longer. Reasons for the different time scales are the two main mechanisms that lead to phase transitioning between gas and brine phase due to chemical potential: dissolution across the pore-scale interfaces in volumes of porous medium containing both gas and brine phase, and dissolution across the macroscopic interface that separates bulks of gas and brine phase, that is, the boundary between gas-phase plume and brine region below. For a representation of the pore-scale and macroscopic interfaces, see Figure 1. Everywhere where gas and brine phase are both present in a volume of porous material, phase transition is an almost instantaneous process, since the pore-scale interfaces separating the phases have a very large specific area. Dissolution across pore-scale interfaces occurs mainly near the injection region, where the gas and brine phases are not in equilibrium and distributed over the entire depth of the aquifer. Furthermore, dissolution across pore-scale interfaces occurs within the developing gas-phase plume, until the residual brine phase disappears due to evaporation of all water into the gas phase (drying front). Dissolution across pore-scale interfaces also plays a role when the mobile gas phase retreats, leaving residual gas phase, which may come in contact with fresh brine, until all gas components are dissolved in brine (wetting front).
Dissolution across macroscopic interfaces is a slow process in general, since the macroscopic interface separating the phases has a relatively small specific area, and the dissolution rate is often limited by diffusion carrying components away from the macroscopic interface. As a result, dissolution across macroscopic interfaces may only be relevant in later times. However, in some cases, it may be accelerated by convective mixing. Convective mixing means that density differences in the phases due to different compositions lead to fluxes within one phase that increases the concentration difference over a macroscopic interface. As an example from CCS, CO 2 that diffuses from the gas-phase plume into the brine phase below increases the brine-phase density. The heavier brine phase sinks to the bottom of the aquifer (often inducing instabilities in form of fingering), bringing fresh, non-saturated brine phase back up to the gas-phase plume. This way, the concentration difference across the boundary of the gas-phase plume is maintained and dissolution of CO 2 remains at a steady rate. This system has been studied extensively in the context of CO 2 storage (  Three different model concepts for compositional VE with increasing complexity: (a) assuming chemical equilibrium over the entire depth of the aquifer, (b) assuming chemical equilibrium only in the gas-phase plume, (c) assuming chemical equilibrium in the gas-phase plume and a rate q of components dissolving into brine below the gas-phase plume.
Note that the height of the gas-phase plume may differ between (a) and (c). The top of the aquifer is denoted by z T and the bottom by z B , while s b represents the saturation of the brine phase.

of 19
Based on the different levels of complexities discussed above, three different model concepts for compositional VE can be identified, as illustrated in Figure 1. Figure 1a shows a simple model concept for compositional VE where the fluid phases and the components are assumed to be in equilibrium over the entire vertical direction. In reality, this state will only be reached in the entire domain after a very long time that is possibly much larger than the time required to reach phase equilibrium. As a result, compositional VE models assuming chemical equilibrium over the entire depth of the aquifer are only applicable for large time scales. If smaller time scales are of interest, compositional VE models may assume chemical equilibrium in the gas-phase plume only, as depicted in Figure 1b. Diffusion/dissolution over the interface between gas-phase plume and brine region below is neglected, assuming pure brine below the gas-phase plume. Compositional VE models assuming chemical equilibrium restricted to the gas-phase plume are of medium complexity since gas-phase plume and region below are treated differently. Applicability of such models is given only for short time scales or cases where there is no convective mixing. Figure 1c shows a model concept of high complexity that deals with the transition between both compositional VE states by assuming chemical equilibrium in the gas-phase plume and applying a component flux across the macroscopic interface. Diffusion/dissolution over the macroscopic interface can be modeled in the form of a component flux that is driven by the concentration difference, or as an upscaled dissolution rate, which may be assumed to be constant. Furthermore, the presence of gas components below the gas-phase plume could also be caused by lateral transport of brine containing gas components or by a gas phase present in this region prior to reaching phase equilibrium. Assumptions have to be made about the profile of the concentration of gas components in the vertical direction below the gas-phase plume, for example, constant concentration of gas components in brine directly below the gas-phase plume, or an averaged concentration over the entire depth below the gas-phase plume.
The following section will transfer the conceptual model of compositional VE to the mathematical model by deriving the governing equations through vertical integration.

Mathematical Models
We first present the three-dimensional governing equations for compositional multiphase flow in a porous medium, formulated as a volume balance that determines the pressure and subsequently solved transport equations for the component concentrations. We define the assumptions of our efficient VE model to enable a reconstruction of the fine-scale solution. Then we derive the coarse-scale and fine-scale equations for the compositional VE model by vertically integrating the three-dimensional governing equations.

Full-Dimensional Model
The component transport equation for each component κ balances the total concentration of the components c κ [M/L 3 ]: where t is the time [T], is the mass fraction [-] of a component κ in phase α, ϱ α is the density [M/L 3 ] of a phase α, u α is the Darcy velocity [L/T] of a phase α, and q κ is the source/sink term [M/(L 3 T)] for each component κ. The phase α can either be the brine phase with a subscript "b" or a gaseous phase with a subscript "g." The extension of Darcy's law for multiphase flow is taken to compute u α : where λ α = k r,α /μ α is the mobility of a phase α We use a volume balance formulation for the pressure equation, which allows to directly include changes of density due to changes of composition. Balancing the fluid volume instead of mass was first proposed by Acs et al. (1985), and similar pressure equations can be found in Watts (1986) and Coats (2000). We assume a BECKER ET AL.
10.1029/2021WR030990 5 of 19 rigid solid matrix for simplicity of notation, neglect diffusion, and refer to Van Odyck et al. (2009) andFritz et al. (2012) for a derivation of the following pressure equation: where is the total specific volume [-] of all fluid phases, p is either one of the phase pressures p α , and ϕ is the porosity [-]. Here, the total specific volume is the sum of specific fluid volume, which describe the relation of the volume occupied by a phase to the total cubature of the control volume. Equation 3 balances volume changes due to compositional changes with changes in the pressure field. The equations of state to determine properties like phase density, phase viscosity and phase equilibrium will be summarized in later parts of this study. The system of equations can then be closed by eliminating one of the two pressures from the equation using the capillary pressure relationship p c = p g − p b , which is assumed to be a function of the brine phase saturation s b [-].
The right-hand side of the pressure equation is called the residual. It is derived from a Taylor expansion of first order of the constraint that the sum of specific phase volumina has to equal the porosity, as presented by Trangenstein and Bell (1989). This residual is introduced because it leads to higher stability in the sequential scheme, where changes in fluid density due to partial miscibility and compressibility are not implicitly incorporated. The time stepping criterion is based on a CFL condition which requires that the representative flux of each cell does not exceed the available storage of that cell. The representative flux is defined as the maximum volumetric outflux or influx of each cell, as proposed by Helmig et al. (2010). All secondary variables in pressure and transport equations are used from the previous time step, and updated at the end of each time step.

VE Model Assumptions
To develop an efficient and robust model for underground energy storage, and in particular to enable reconstruction of the fine-scale solution in the vertical direction, a number of assumptions and simplifications for our compositional VE model are necessary. The components of interest in underground energy storage (e.g., hydrogen and methane) do not induce a convective mixing behavior in our simulations. This is because these components have smaller molar masses than H 2 O and would decrease the density of the brine, according to simplified mixing laws that assume that one molecule of liquid component is replaced by one molecule of gas component (Class & Helmig, 2002). Consequently, our compositional VE model neglects the dissolution of gas components below the gas-phase plume and assumes that only the region of the gas-phase plume is in chemical equilibrium (VE state (b) in Figure 1). We note that, if necessary, the model can be extended to include dissolution across the macroscopic interface separating gas-phase plume and brine phase below in the same way as presented by Gasda et al. (2011b), Nilsen et al. (2016), or Møyner et al. (2020. We largely neglect the influence of pressure variations along the vertical direction on equilibrium mass fractions, phase densities, and phase viscosity, by assuming a constant value for these properties inside the gas-phase plume and below, depending on reference pressures inside the gas-phase plume and below. The reference pressure varies in location among the properties to be reconstructed, but is always a fine-scale pressure, which is a reconstructed pressure on the full-dimensional discretization. For gas-phase properties, the reference pressure is the gas phase pressure at the top of the aquifer. For properties of the brine phase inside the plume, the reference pressure is the brine phase pressure at the lower end of the gas phase plume while for properties of the brine phase below the plume, the reference pressure is the brine phase pressure at the bottom of the aquifer. These simplifications should not introduce significant errors since property variations along the vertical direction are less significant than property variations along the horizontal direction, as long as the vertical extent of the domain is much smaller than the horizontal extent (Bachu, 2002). Furthermore, the variation of gas-phase pressure along the vertical direction is expected to be small due to the low density of the gas and the brine phase density will be close to constant in the vertical direction due to the brine phase being almost incompressible. In addition, we assume that the plume region is a two-phase region in case residual saturations are larger than zero (or for cases with a capillary transition zone) and consequently neglect drying and wetting fronts. 6 of 19

Reconstruction
With the model assumptions above, we can reconstruct the fine-scale solution in the vertical direction, as shown qualitatively in Figure 2 for the brine phase saturation, phase pressures, mass fractions of the water component, and phase densities. The profile of the brine phase saturation in the gas-phase plume follows the inverse of the capillary pressure curve and is piecewise constant below the gas-phase plume. The pressure distribution is assumed to be hydrostatic, with different constant slopes in the region of the gas-phase plume and below it. Densities and mass fractions (and, analogously, viscosities, which are not pictured) inside the gas-phase plume and below are piecewise constant, due to different compositions and different reference pressures in these regions. We note that during reconstruction of the vertical saturation profile, we fill up the cell from the top with the volume of both phases to determine the height of the gas-phase plume. The reconstruction mechanism hereby ensures that gas components and associated volume is consistent with the volume of the available pore space. As a result, the volume mismatches allowed by the compositional model can only be caused by water components. This neglects the volume mismatch in the plume region of a VE cell and shifts it to the bottom of the VE cell. A volume mismatch is the discrepancy between pore volume and fluid volume, which is caused by the sequential solution scheme and is reflected by the residual term in Equation 3.

Governing Equations
The governing equations of the compositional VE model are derived by integrating the full-dimensional equations over the height of the aquifer from the bottom of the aquifer, z B [L], to the top of the aquifer, z T [L]. We identify the reconstructed, detailed solution in the vertical direction as the fine scale and the vertically integrated equations as the coarse scale (Nordbotten & Celia, 2011). We refer to vertically integrated variables as coarsescale variables, denoted by uppercase letters, and to the variations along the vertical direction as fine-scale variables that are denoted by lowercase letters. First, we integrate the pressure Equation 3: Due to the assumption of hydrostatic pressure distribution, is constant over z and neglecting the influence of pressure, is piecewise constant inside the gas-phase plume and below it, respectively. Further, we move the volume derivatives ̂ and ̂ out of the integral because we obtain these terms as a representative constant value for the entire VE cell, by varying P and C κ and observing the change of ̂ after reconstruction. This leads to: with ̂= ∫̂ , = ∫ , Φ = ∫ and = ∫ . The subscript "//" denotes the xy-plane, the superscript "VE" indicates that we obtain the values in parenthesis during VE reconstruction. Properties below the gas-phase plume are denoted with superscript ( * ), inside the gas-phase plume with superscript ( * * ). The coarse scale Darcy velocities * and * * are defined as: with the minimum vertical location of the gas-phase plume min [L], that stores the minimum value of z p (bottom of mobile gas-phase plume) over time. * is defined at the bottom of the aquifer, while * * is defined at the top of the aquifer. The depth-integrated permeability and depth-averaged relative permeability below and above the gas-phase plume are defined as: The VE transport equations are determined by integrating Equation 1, as for each component κ. The pressures * * at the top of the aquifer are related to the pressures * at the bottom of the aquifer by exploiting the assumption of hydrostatic profiles and the minimum location of the gas-phase plume known from the earlier time level. To solve the equation, one of the two remaining phase pressures at the bottom of the aquifer is then eliminated using the (pseudo) capillary pressure = * − * .

Hybrid Framework: Coupling VE and Full-Dimensional Model
The compositional VE assumption may not be valid everywhere in the domain, mainly due to the phases not having reached equilibrium yet. This especially concerns the vicinity of injection areas, where injected gas forms a gas phase that moves upward, and the tip of the advancing/retreating gas-phase plume, where brine drains out of or imbibes the gas region. In proximity of extraction areas, the detailed location and composition of the gas phase are of high importance, which is why a full-dimensional model may be desired here, too. Therefore, we couple the compositional full-dimensional model and the compositional VE model in one domain. In the following, we first explain the coupling by formulating the volume change due to fluxes across the subdomain interfaces in Section 4.1, and present the computational algorithm to solve the governing equations of the hybrid model in Section 4.2.

Volume Change Due To Fluxes Across Subdomain Interfaces
We exploit the similarity between the pressure Equation 3, including Equation 2, of the compositional full-dimensional model, and the pressure Equation 5, including Equations 6 and 7, of the compositional VE model to couple both models in one domain. Both equations balance a storage term describing the change of total specific volume of all phases due to pressure changing over time, with volume changes due to fluxes, a source/sink term and a residual describing discrepancies in total specific volume and porosity. The two compositional models are coupled together by formulating the change of total specific volume of phases due to fluxes across the subdomain interfaces (second term in Equations 3 and 5) and requiring continuity. For that purpose, we refine the VE cell on one side of the subdomain interface into full-dimensional subcells in the vertical direction (see gray dots and dotted lines in Figure 3). Each subcell shares an interface with a neighboring cell in the full-dimensional subdomain.
The discretized form of the volume change due to fluxes from the full-dimensional cell i to a VE subcell j′ on the other side of the boundary between subdomains is: where all coefficients without the superscript t are calculated on the previous time step, whereas t indicates the current time step. The normalized vector n γ is orthogonal to the interface γ shared by cells i and j′ and pointing outward of the cell i, d ij′ is the normalized vector connecting the centers of the cells i and j′ and pointing away from i, and ′ is the face area shared by the full-dimensional cell and the VE subcell. The permeability tensor k′ is determined by calculating the harmonic mean of corresponding entries of the tensors in the full-dimensional cell and the VE subcell, while V i describes the cell volume and U i the cell surface. The quotient in the last summation over the components κ in Equation 13 is based on the assumption that the gradient of each volume derivative can be approximated linearly between neighboring cells. However, the term's derivation is not straightforward and lengthy, therefore we refer to Fritz et al. (2012). One phase pressure p α,j′ , in our case p g,j′ , is replaced by the sum of the other phase pressure and (̄′ ) , where ′ is the brine phase saturation in the VE subcell that would form under conditions of chemical equilibrium within the subcell, using averaged total concentrations and reconstructed pressures of the subcell from the previous time step to establish the state of chemical equilibrium. The remaining phase pressure p α,j′ is the reconstructed pressure at the calculation point of the VE subcell and is expressed using the phase pressure at the bottom of the VE cell, * . In general, we account for a capillary transition zone inside the gas region, so brine can flow within the transition zone. Furthermore, we use full upwinding for the reconstructed flux-related properties phase density ′ , phase mobility ′ and mass fraction ′ , depending on the phase potential. The phase potentials are: where ′ is centrally averaged between cell i and j′, and pressures and densities are taken from the previous time step. If the sign of the potential is positive, the flux-related properties for phase α are taken from cell i; if it is negative, from the VE subcell j′. All flux-related properties in the VE subcell j′ are calculated assuming chemical equilibrium within the subcell, using averaged total concentrations and reconstructed pressures of the subcell from the previous time step. Furthermore, the volume derivative ̂ ′ is centrally weighted between cells i and j′, while ̂ ′ is determined specifically for the subcell itself, using averaged properties inside the subcell and a fraction of the transport estimate according to the volume ratio between subcell and VE cell. For each interface of a full-dimensional cell, there is only one adjacent VE subcell contributing to volume change due to fluxes from the VE subdomain. The volume change caused by fluxes over the subdomain interface to the VE cell is computed as the sum of the volume changes due to fluxes from all neighboring full-dimensional cells, while for a full-dimensional cell only the contribution of the adjacent subcells is taken into account.

Computational Algorithm
We solve the pressure and transport equations of the hybrid model sequentially, using an IMPET scheme. Before each time step calculation, we compute a transport estimate on the pressure field of the previous time step to determine the volume derivatives. Then, we first solve the pressure equation in an implicit manner, where the primary unknown variable is one phase pressure for each computational cell, and all coefficients are calculated based on the previous time level. We use a single computational matrix for the entire domain, combining full-dimensional subdomains (pressure Equation 3, including Equation 2) and VE subdomains (pressure Equation 5, including Equations 6 and 7) without iterating between subdomains, and consequently solve for all phase pressures simultaneously. The volume change due to mass fluxes at the interface is calculated as described above, with all adjacent full-dimensional cells taken into account for a VE cell at the subdomain interface.
Once the pressure has been updated, the full-dimensional transport (Equation 1) and VE transport (Equation 12) are solved explicitly to update total concentrations in each cell. The scheme is mass conservative since all mass leaving one cell enters the neighboring cell. Once the total concentrations are known, all secondary variables are updated.

Adaptivity
During the simulation, areas in the domain that are not in VE may change in size and position, depending on local processes like advancing gas-phase plume or geological heterogeneity, or because phase equilibrium is only reached after a certain time. Therefore, during simulation, our model identifies the parts of the domain that are not in equilibrium, and adaptively assigns the compositional full-dimensional model in these parts. Everywhere else, the compositional VE model is used. For that purpose, we develop an a posteriori criterion for local VE in the context of compositional models in Section 5.1, and analyze the behavior of the criterion value using a test case of gas storage in Section 5.2. Finally, Section 5.3 presents the adaptive algorithm that assigns subdomains to the compositional VE and full-dimensional model.

Local VE Criterion
We develop a local a posteriori criterion for VE in the context of compositional models. We chose a criterion based on profiles of relative permeability in the vertical direction, as proposed by Becker et al. (2018) for immiscible two-phase flow. In this previous work, the criterion values showed promising behavior over space and time with a distinct peak at the leading edge of the gas-phase plume, which simplifies setting a threshold value for switching between models. We note that pressure profiles could also be used for the criteria, though we have not explored this in much depth. Furthermore, relative permeability profiles govern overall flow behavior, since the fine scale relative permeability is integrated to receive the coarse scale relative permeability. Thus, relative permeability profiles are a natural choice to determine applicability of VE models. To receive the criterion, the difference between the simulated profile in a full-dimensional column to a hypothetical VE profile is calculated as the area between the profiles, and normalized with the thickness of the mobile gas-phase plume in the hypothetical VE profile (Becker et al., 2018). The equation to calculate the criterion values thus follows as: where ′ is the reconstructed relative permeability and h is the thickness of the mobile gas-phase plume in the hypothetical VE profile with h = z T − z p . The criterion computed in Equation 15 equals the shaded area in Figure 4 scaled with h. Relative permeability is reconstructed by applying a relationship like, for example, Brooks-Corey, to a reconstructed profile of brine saturation. Brine saturation is reconstructed assuming that all gas components inside the full-dimensional grid columns are located in the two-phase region of the gas-phase plume and that the phases in this region are in chemical equilibrium. Phase densities and mass fractions are estimated according to the phase pressures in the full-dimensional column at the locations where the reference pressures are evaluated. As a result, the values of the developed criterion are directly associated with phase equilibrium, while chemical equilibrium within the gas-phase plume is indirectly included through the reconstruction, to the degree that chemical equilibrium is relevant for the overall flow behavior. Thus, the criterion values can be used to decide if phase equilibrium or the assumption of chemical equilibrium in the gas-phase plume is violated, for example, by a large amount of gas components dissolved below the gas-phase plume. If the criterion does not hold anymore, VE columns can be separated into stacks of cells and vice versa.
We extend the local criterion to take gas extraction, as well as injection, into account, as required for gas storage scenarios. As the gas-phase plume retreats during gas extraction, it leaves residual gas phase, which is in chemical equilibrium with the surrounding brine, or, in case of zero residual gas saturation, a region of brine with dissolved gas at equilibrium conditions. The minimum location of the gas-phase plume, min , is a parameter that depends on the system's history, which can be viewed as a hysteretic behavior on the coarse scale (Doster et al., 2013a(Doster et al., , 2013bJuanes et al., 2010). The mobile gas phase above the immobile region may be in equilibrium, which should be reflected by the VE criterion. For this purpose, we store min for each full-dimensional grid column, as calculated during the reconstruction of the hypothetical VE saturation profile. When extraction occurs, the hypothetical VE saturation profile takes the stored min into account, assigning residual saturation and chemical equilibrium below the mobile gas region. Figure 4 shows an example of simulated and hypothetical VE saturation profiles in a full-dimensional column during extraction, as well as the simulated and hypothetical relative permeability profiles that are used to calculate the VE criterion.

Criterion Analysis
We analyze the behavior of the VE criterion values over space and time for an injection and an injection/extraction scenario with CH 4 in 1,000 m depth, as depicted in Figure 5. Gas is injected from the left-hand side over the entire thickness (30 m) of an initially brine-saturated domain. In the injection/extraction case, gas is then extracted from the left-hand side over a vertical section of 5 m below the top of the aquifer. The domain is long enough such that the gas components will not reach the right-hand side boundary during the . The blue curve is the result of the full-dimensional solution, and the orange curve is the reconstructed profile that would develop if the fluids were in VE. The difference between the relative permeability profiles is depicted using the striped areas. As these profiles belong to a column close to an injection/extraction area at early times, the criterion probably is not fulfilled in the case above. simulation. For the injection case, we uniformly inject 0.0175 kg s −1 m −1 CH 4 for the duration of the simulation (240 h). For the injection/extraction case, we use the same injection rate for 96 h and then we uniformly extract methane at a constant extraction rate of 0.003 kg s −1 m −1 . We use Brooks-Corey curves for relative permeability and capillary pressure with pore size distribution index λ = 2 and entry pressure p e = 1.0 × 10 5 Pa. The residual saturation of the brine phase is zero while the residual saturation of the gas phase is assumed to be 0.1. We use a grid resolution of 1.0 m in horizontal and 1.0 m in vertical direction. We calculate the segregation time after Nordbotten and Dahle (2011).
where H is the height of the aquifer [L], s b,r is the residual brine phase saturation [-], and k z represents the vertical permeability [L 2 ]. The segregation time for our scenario is estimated as t seg = 48 h. In theory, the relative permeability k r,b should be a value which is a representative for the entire process of segregation. Since that is essentially unknown, in practice it is usually set to k r,b = 1.0.
In our compositional simulations with CH 4 , we use the following equations of state to determine the phase density and phase viscosity, as well as the Henry coefficient and vapor pressure for the phase equilibrium. We assume local chemical equilibrium in the computational cells of the full-dimensional model and apply a fugacity-based approach assuming ideal mixtures and an ideal gas phase with regard to mixing. Gas and brine phase densities are calculated depending on pressure, temperature, and phase composition. Density of the pure gas phase is calculated assuming the ideal gas law, density of the pure brine phase is calculated according to the IAPWS formulations (IAPWS, 1997). The density of the gaseous mixture is calculated assuming an ideal mixture, while the density of the liquid mixture is calculated assuming that one molecule of water component is replaced by one molecule of gas component (Class & Helmig, 2002). The Henry coefficient for the gas component in brine is calculated depending on temperature according to the IAPWS formulations with coefficients published in Watanabe and Dooley (2004). The vapor pressure of pure brine is calculated according to the IAPWS formulations as well, depending on temperature. Phase viscosities depend on temperature and pressure, as well as on composition for the gas phase. Viscosity of the pure gas phase is calculated according to the formulation by Reid et al. (1987), using coefficients published by Poling et al. (2001). The viscosity of the gaseous mixture is calculated using the Wilke method from Reid et al. (1987). Viscosity of the brine phase is calculated according to the IAPWS formulation for water (Cooper & Dooley, 2008).
For the injection case, Figures 6a and 6b show snapshots of the values of the relative permeability criterion, as calculated in Equation 15, in the whole domain for three times, and the behavior of the criterion values over dimensionless time for two different locations along the horizontal axis. Similarly, as observed for two incompressible, immiscible phases by Becker et al. (2018), the criterion values show a high peak around the injection region, indicating that the phases here are not in equilibrium due to injected gas-phase moving continuously upward. Toward the tip of the gas-phase plume, the criterion values increase again, indicating that brine phase drains out of the advancing gas-phase plume. Over dimensionless time the criterion values show a nonmonotone behavior, again similar to the case of two incompressible, immiscible phases. This is due to the finite size of the cells in the full-dimensional model and loses impact with time and a growing thickness of the gas-phase plume as a result.
For the injection/extraction case, Figures 6c and 6d show that over the length of the domain, extraction leads to a slight increase of the criterion values in a region around the tip of the gas-phase plume. This is due to the thickness of the gas-phase plume decreasing, especially toward the tip, causing the criterion to become more sensitive, as it is normalized by the thickness. At a fixed location inside the gas-phase plume, the criterion values stay constant over time as soon as extraction starts. This is opposed to the time before extraction and the case without extraction, where the criterion values decrease with time. During extraction, the effect of decreasing thickness of mobile gas-phase plume and the system striving toward VE cancel each other out. Close to the injection/extraction area, where most of the gas is extracted, the criterion values even increase at later times, indicating that here the influence of small gas-phase saturations below the gas-phase plume is dominating.
We note that our compositional VE criterion values may show a decrease over time for the injection case, indicating that VE is approached, but we expect the criterion values to increase again after very long time scales. This is the time scale after which the system moves toward chemical equilibrium along the vertical height of the aquifer, and a compositional VE model assuming chemical equilibrium only in the gas-phase plume becomes unsuitable. In this case, gas components dissolving across the macroscopic interface between gas-phase plume and brine phase below should be included in the model. This did not become critical in our simulations, mostly because the process of dissolution across this macroscopic interface is not accelerated by convective mixing in our full-dimensional model and the effect on the VE criterion is counteracted by the phases striving toward phase equilibrium. Thus, the assumption of chemical equilibrium in the gas-phase plume remains sufficiently valid throughout our simulation times.
Choosing a threshold value for the VE criterion ϵ crit is not always straightforward. On the one hand, it has to be chosen manually, and on the other hand, it needs to be an efficient value. For a comparison of different threshold values and their impact on results, see Becker et al. (2018). Especially, a scenario with several injection/extraction cycles for a time long enough for most of the plume to be in VE poses a challenge. Then, it may just happen that the VE region around the tip vanishes during one of the injection cycles. It never returns afterward. At first sight, that does not seem as a problem since the entire plume is in equilibrium. However, if the injection rate or time is increased and the tip moves into a region where the plume has not been before, while reaching a lens or some other vertical structure requiring full-dimensional simulations, it will completely miss it. All in all, losing the tip when conditions change after that, is a problem which is ultimately based on the initial choice of the criterion threshold.

Adaptive Algorithm
The decision as to which parts of the domain are in VE is made based on the VE criterion, and is evaluated after each time step for all full-dimensional grid columns. A full-dimensional grid column where the VE criterion value is smaller than a threshold ϵ crit is considered to be in VE. The following decision is made for each grid column: 1. A full-dimensional grid column stays full-dimensional if the VE criterion is not met (c relPerm ≥ ϵ crit ). 2. A full-dimensional grid column is turned into a VE column if the VE criterion is met (c relPerm < ϵ crit ), and the column is not a direct neighbor to a column where the criterion is not met.
BECKER ET AL.
10.1029/2021WR030990 13 of 19 3. A VE cell is turned into a full-dimensional grid column if it is a direct neighbor to a (full-dimensional) column where the criterion is not met.
The second and third requirements lead to a buffer zone around regions that are not in VE. The buffer cells serve as indicators to turn VE cells back into full-dimensional grid columns, as done by the third requirement. Due to the CFL restriction of the sequential solver, the solution does not propagate more than one spatial cell within one time step, suggesting that the buffer zone should contain one layer. However, the buffer zone can be extended to have more than one layer for stability reasons. In addition, we recommend that information that is known beforehand, for example, injection at a later point in time in another region of the domain, or a fault zone where vertical movement of brine is of interest, be implemented by requiring full-dimensional cells in this region.

Results and Discussion
We test the applicability of the compositional VE model and the compositional hybrid model using a test case of gas storage in an idealized domeshaped aquifer with two low-permeability lenses. The test case includes primary injection of gas into a previously brine-saturated heterogeneous aquifer, and long-term development of the gas-phase plume after several injection and extraction cycles. The compositional hybrid model, the compositional VE model, and the compositional full-dimensional model are all implemented in DuMu x (Flemisch et al., 2011;Koch et al., 2020), allowing for a comparison of the different methods within the same software framework.

Gas Storage in an Idealized and Heterogeneous Underground Aquifer
We simulate hydrogen storage in a radially symmetric anticline structure consisting of sandstone in 1,700 m depth, see Figure 7. Additionally, two lenses with a low permeability (0.022mD) are added, located at 30 and 120 m distance from the injection location, respectively, with a length of 10 and 20 m, respectively, and with a thickness of 10 m each. The aquifer is discretized into a radially symmetric grid with horizontal resolution of 5 m and vertical resolution of 0.94 m (for full-dimensional model and full-dimensional subdomains of the hybrid model). A threshold value of 0.07 for the relative permeability criterion was chosen, starting with a higher value which is reduced until a stable solution is reached, as suggested by Becker et al. (2018). Hydrogen is injected and extracted in the center of the radially symmetric aquifer, through a large-diameter well (0.53 m) that is perforated in the top 5 m of the aquifer. We use Brooks-Corey curves for relative permeability and capillary pressure with pore size distribution index λ = 2 and entry pressure p e = 1 × 10 5 Pa. The residual saturation of the brine is zero while the residual saturation of the gas phase is assumed to be 0.1.
The operation of the subsurface storage is split into two parts. During a 30-month development period, we inject 1.571 × 10 −3 kg s −1 H 2 for 3 months, followed by seven idle months, repeated three times. During the operation period, 7.845 × 10 −4 kg s −1 of H 2 are extracted for 3 months, followed by injection of 7.845 × 10 −4 kg s −1 of H 2 for 3 months, for a total of 24 months. The segregation time is estimated as t seg = 10 months. Figure 8 shows the injection rate, the total mass of injected gas since the start of the simulation, and the bottom hole pressure, for the full-dimensional simulation. The pressure response during the development period shows that because of the open boundary on the far side of the well, the pressure buildup due to injection is lost within approximately 1 month of idle time. Figure 9 shows the total concentration of gas components for all three models after the first injection during the operation period, as well as after the last extraction and the last injection during the operation period. The hybrid and the full-dimensional model are in very good agreement considering plume shape, especially around the heterogeneities and regarding the area of residual gas phase saturation, throughout the simulation. The VE model assumes segregated phases from the beginning, which results in a very different plume shape in early times, and Figure 7. Cross section of the geometry for the hydrogen storage structure with initial and boundary conditions as well as the material properties. To determine density and viscosity of hydrogen as well as phase equilibrium, we use the Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) (Gross & Sadowski, 2002). Sauerborn (2021) showed that appropriate treatment of density and viscosity of hydrogen has an influence on pressure development close to the injection point of a hydrogen storage aquifer. For water, the highly specialized and verified IAPWS formulations (IAPWS, 1997)   completely misses the low-permeability lenses in later times, as depicted in the top row of Figure 9. The general plume length is predicted similarly by all models. We note that for the VE model, this is due to a combination of the low-permeability lenses filling with gas, which slows down the advancing plume, and assuming the phases have segregated, which accelerates the advancing plume, so that in total the plume length is only sightly overestimated. The full VE model predicts a different plume shape in the vicinity of the well than the other two models, with a larger plume height after injection and a smaller plume height after extraction directly at the well. The total gas concentration in the gas plume is slightly higher for the full VE model and the VE subdomains of the hybrid model, while the thickness of the gas plume is smaller compared to the full-dimensional model. This is due to the finite grid discretization of the full-dimensional model, which leads to a slight smearing of the phases over the vertical direction. Figure 10 compares pressures in the computational cell directly at the bottom of the well, composition and saturation of gas phase in the computational cell directly at the bottom of the well, and total mass of dissolved gas components in the domain for all three models. In terms of pressure and composition (Figures 10a and 10b), the hybrid model compares very well to the full-dimensional model close to the well. The VE model slightly underestimates the pressure increase in the domain close to the well, due to assuming hydrostatic pressure distribution in the vertical direction, which is violated in the area close to the well. As a direct effect of pressure deviations at that location, the VE model slightly overestimates the mass fraction of water components in the gas phase during injection. The gas-phase saturation close to the well (Figure 10c) is captured equally well by the full-dimensional and the hybrid model, while the VE model shows strong deviations because the gas and brine phase are not in equilibrium near the well. Figure 10d shows the total mass of dissolved gas components in the domain over time for the three models. The full-dimensional model and the adaptive hybrid model show a steep gradient in mass of dissolved gas components

of 19
at the beginning of the simulation. This strong increase in dissolved gas is caused by injection of gas into a fully brine-saturated aquifer, where it dissolves immediately until equilibrium concentration is reached. Only then a separate gas phase forms, and moves upward due to buoyancy. This effect can be observed even though we only inject over the top 5 m of the aquifer instead of the entire vertical length. This process is not captured by the full VE model, which assumes from the beginning that most gas has formed a separate phase and moved to pool below the top of the aquifer. The overall trend of the full-dimensional model shows an increase of the mass of total dissolved gas components during the operation period, while the hybrid model and in particular the VE model predict no overall increase during the operation period. At later times, both the VE model and the hybrid model underestimate the total mass of dissolved gas components predicted by the full-dimensional model. However, due to hydrogen's low solubility, the difference in total mass of dissolved gas components computed by the full-dimensional model and the hybrid model at the end of the simulation only makes up 0.6% of the total mass of gas components in the aquifer.
In general, the VE model and VE subdomains underestimate the mass of dissolved gas due to assuming one reference pressure in the gas-phase plume to approximate equilibrium conditions. This reference pressure is chosen as the gas phase pressure at the top of the aquifer, that is, the smallest gas phase pressure along the vertical direction. Moreover, the VE model and VE subdomains neglect dissolution of gas components into the pure brine below the gas-phase plume, whereas the full-dimensional model develops a small region with gas-saturated brine around the gas-phase plume. For longer simulation times than studied here, dissolution of gas components into the brine below the gas-phase plume could be included in the model by switching to fully gas-saturated brine below the gas-plume in VE subdomains (model in Figure 1a).
Finite discretization in the vertical direction has an influence on the total mass of dissolved gas components, as well. Simulation results of a full-dimensional model on a vertically refined grid (vertical resolution of 0.23 m) are closer to the hybrid model, as shown in Figure 10d, dash-dotted line. Finite discretization in the vertical direction causes a smearing of the gas phase along the vertical direction, which leads to more gas components dissolving into brine. Additionally, for the full-dimensional model, we observe very small concentrations of dissolved gas that are transported away from the gas-phase plume during recurring injection and extraction phases. This effect reduces with refinement of the grid.

Performance Comparison
On average, during the simulation, the hybrid model uses 42% of the computational cells required for the full-dimensional simulation and 36% of the CPU time compared to the full-dimensional model. The full VE model uses 0.4% of the CPU time. Since the efficiency of the hybrid model is directly related to the ratio of VE subdomain area to full-dimensional subdomain area, we expect the hybrid model to become even more favorable in the case of real scenarios, where large parts of the domain are either one-phase regions or occupied by the gas-phase plume in VE.

Influence of Lens Permeability
We increase the permeability of the lenses in the above storage case from the original 0.022-0.22 and 2.2mD, respectively, and analyze the influence of lens permeability on the hybrid model. Figure 11 shows the total concentration of gas components at the end of simulation for the hybrid model with a relative permeability criterion of 0.07 and for the full-dimensional model, for all three lens permeabilities. In the case with the highest lens permeability, only the injection region and the tip of the gas-phase plume are assigned full-dimensional subdomains at the end of the simulation, since both lenses have filled with gas-phase resulting in phase equilibrium in the vertical direction. For the case with a medium lens permeability, the smaller lens is filled with gas phase at the end of the simulation, while the larger lens is still not in VE, according to the relative permeability criterion. Consequently, the injection region, the tip of the gas-phase plume, and the larger lens are full-dimensional subdomain at the end of simulation. In comparison with the full-dimensional model, it becomes clear that the region of residual gas-phase saturation around the injection region and around the smaller lens is not represented well in those cases where the smaller lens is turned into a VE subdomain, see the second and third column of Figure 9. Here the model fails to capture effects from the first injection periods, where the gas phase was not in equilibrium with the brine phase and spread over a larger vertical extent, causing a larger area of residual saturation than predicted by the VE model. When the model switches from full-dimensional to VE subdomain, these extended areas of residual gas phase saturation are lost. As in all cases where the point of injection is known prior to simulation, we recommend to manually prescribe a larger full-dimensional subdomain around the injection region, where such effects are to be expected. We note that to some extent, the smearing of the gas-phase across the vertical direction and the resulting region of residual gas phase may also be caused by the finite discretization in the vertical direction.
Independent of the lens permeability, all simulations start out with the entire plume as full-dimensional subdomain in early times. As the plume advances, both lenses are assigned full-dimensional subdomains at early times in all simulations as well. At later times and depending on the permeability of the lenses, VE subdomains appear inside the plume at different times and in different sizes. Figure 12 shows the number of computational cells over simulated time for all three lens permeabilities. For the case with the highest lens permeability, one can clearly observe the drop in the number of computational cells after approximately 10 months, when the smaller lens near the injection region turns into a VE subdomain. Next, after approximately 3 years simulated time, the larger lens turns into a VE subdomain as well, causing another decline in the overall number of computational cells. For the case with the medium lens permeability, the smaller lens near the injection region turns much later into a VE subdomain, after around 3 years simulated time. This shows that the hybrid model adapts very well to the different local segregation times, and that with time the simulation becomes more efficient as more and more regions of the domain turn into VE subdomains as the phases segregate.

Conclusion
We have developed a hybrid VE/full-dimensional modeling framework for compositional multiphase flow in the context of subsurface energy storage, which dynamically determines and couples VE and full-dimensional domains. For this purpose, we derived a compositional VE model by vertical integration of the transport equation and volume balance equation of the compositional full-dimensional model. The compositional VE model showed good accuracy in predicting the long-term gas-plume length compared to a full-dimensional model, but missed local heterogeneities and details in the gas-plume shape. Furthermore, it showed low accuracy in predicting the gas-phase saturation close to the injection area. The compositional hybrid model showed excellent accuracy in predicting the gas plume development, as well as composition and pressure development near the well with a lower number of computational cells and consequently lower computational cost compared to a full-dimensional model. The regions that are in VE increase with time as phases segregate, which makes the compositional hybrid model favorable for longer simulation times. The improvement of computational efficiency is a function of the ratio between VE and full-dimensional subdomains -the greater the VE subdomain, the more efficient the hybrid model will be. For the specific test scenario of gas storage in an idealized 2D dome-shaped aquifer used in the present study, the simulations using the hybrid model took 36% of the CPU time compared to that of the full-dimensional model. For larger simulation scales, we assume that the efficiency of the hybrid model will increase even more, as the ratio of VE subregions to full-dimensional subregions should also increase. If highly miscible gases and long simulation times are of interest, dissolution below the gas-phase plume could be included in the VE model and VE subdomains, by assuming fully gas-saturated brine below the gas-phase plume. As such, the compositional hybrid model is an efficient tool for modeling underground energy storage with a focus on local gas-phase saturation and composition.