Numerical Modeling of Gas Hydrate Recycling in Complex Media: Implications for Gas Migration Through Strongly Anisotropic Layers

Burial driven recycling is an important process in the natural gas hydrate (GH) systems worldwide, characterized by complex multiphysics interactions like gas migration through an evolving gas hydrate stability zone (GHSZ), competing gas‐water‐hydrate (i.e., fluid‐fluid‐solid) phase transitions, locally appearing and disappearing phases, and evolving sediment properties (e.g., permeability, reaction surface area, and capillary entry pressure). Such a recycling process is typically studied in homogeneous or layered sediments. However, there is mounting evidence that structural heterogeneity and anisotropy linked to normal and inclined fault systems or anomalous sediment layers have a strong impact on the GH dynamics. Here, we consider the impacts of such a structurally complex media on the recycling process. To capture the properties of the anomalous layers accurately, we introduce a fully mass conservative, high‐order, discontinuous Galerkin (DG) finite element based numerical scheme. Moreover, to handle the rapidly switching thermodynamic phase states robustly, we cast the problem of phase transitions as a set of variational inequalities, and combine our DG discretization scheme with a semi‐smooth Newton solver. Here, we present our new simulator, and demonstrate using synthetic geological scenarios, (a) how the presence of an anomalous high‐permeability layer, like a fracture or brecciated sediment, can alter the recycling process through flow‐localization, and more importantly, (b) how an incorrect or incomplete approximation of the properties of such a layer can lead to large errors in the overall prediction of the recycling process.

Given the vast complexity of the gas hydrate systems, the modeling of gas hydrate dynamics in general, and burial-driven recycling in particular, poses multiple conceptual and computational challenges. One interesting challenge that is central to gas hydrate dynamics is that of gas-water-hydrate phase transitions. In such situations methane can dissolve into and exsolve from porewater leading to a locally appearing and disappearing free-gas phase, and the gas hydrates can melt or precipitate, leading to an appearing, disappearing, and evolving solid phase. These fluid-fluid and fluid-solid interactions occur at different time scales, for example, methane dissolution-exsolution is a fast process governed by vapor-liquid equilibrium on the geological time scales, while gas hydrate phase change is a slower process where the equilibrium assumption may not hold under rapid sedimentation, and is therefore, modeled as a kinetically controlled process (Gupta et al., 2020). These phase transitions are in permanent competition that drives the aforementioned cyclic rebuilding of gas hydrates.
The numerical challenges related to the phase transitions are discussed in (Class et al., 2006;Marchand et al., 2013). Different numerical techniques have been constructed to overcome the phase transitions in multi-phase multi-components porous media models, for example, primary variable switching (PVS) (Class et al., 2002;Wu & Forsyth, 2001), negative saturations (Panfilov & Panfilova, 2014), method of persistent variables (Huang et al., 2015;Neumann et al., 2013) and non-linear complementary problem approaches (Ben Gharbia & Jaffré, 2014;Kräutle, 2011;Lauser et al., 2011). PVS schemes are implemented in many of the hydrate reservoir simulators such as TOUGH-Hydrate (Moridis et al., 2008). Gupta et al. (Gupta et al., 2020) extended the non-linear complementary constraints approach of (Lauser et al., 2011) to gas hydrate systems, and showed that under rapidly switching phase states, this approach seems to be capable of handling the gas-water-hydrate phase transitions more accurately, robustly, and efficiently compared to the more traditional PVS schemes. Furthermore, with the help of this newly developed simulator, Schmidt et al. (2022) were able to demonstrate the mechanics of the hydrate nozzle and its implications for gas migration through GHSZ during continuous burial.
Conventionally, gas hydrate recycling has been largely studied in 1D geological settings and the underlying sediments are assumed to be either homogeneous (Schmidt et al., 2022) or with vertically stacked topography (You et al., 2019(You et al., , 2021) that is representative of the different granular materials, debris, and organic matter that was deposited over different geological times in the past. However, complex fault systems, fluid escape structures, and anomalous sediment layers have been observed in the seismic profiles cross-cutting the buried layers within the GHSZ worldwide (Crutchley et al., 2021;Paganoni et al., 2018;Portnov et al., 2019;Waage et al., 2019). In fact, the formation and propagation of focused flow pathways (pipes, chimneys) and their implications on gas migration, hydrate dynamics, and slope stability remains an important open question. Numerically, handling such anomalous layers is quite challenging and computationally expensive in terms of how the heterogeneity and anisotropy is physically modeled, how the numerical scheme approximates the related fluxes and any material interfaces and/or discontinuities. In the past, we have used the cell centered finite volume (FV) methods for numerical discretization of our gas hydrate models (Gupta et al., 2015(Gupta et al., , 2020 because of their inherent simplicity for implementation of non-linear complementary problems (NCP), local mass conservation property, monotonicity, and low computational costs due to their low order and small two-point stencils. While FV methods offer a very robust, efficient, and reliable numerical framework for simple geological media, it is notoriously difficult to extend to unstructured meshes, fully anisotropic media, and discontinuous material interfaces. Forms of finite element (FE) (Cheng et al., 2013;Fang, 2010) and finite difference (FD) (Holder & Angert, 1982;Yu et al., 2017) methods are also commonly used for methane hydrate models, but they also face challenges related to phase transitions, local mass conservation, overshoots and undershoots (which further complicate the phase change problem), mesh sizes and local mesh anisotropy, and material interfaces. The discontinuous Galerkin (DG) finite element method generalizes the FE method by omitting continuity constraints, allowing potential jumps through numerical fluxes (Cockburn et al., 2000). Moreover, DG methods are locally conservative and a consistent flux across the element interfaces can be easily constructed. Therefore, DG methods, which are generalization of both FV and FE methods, appear to be more suitable for the numerical solution of the methane hydrate model not only because it can handle complex geometries and meshes (including hanging nodes), full material anisotropies, and jumps across material interfaces in a natural manner without additional computational overheads, but also because it preserves the local mass conservation property of the FV method while at the same time provides higher order approximations like FE methods. Moreover, DG schemes offer massive parallelization capability (Bey et al., 1996), which is very important for practical applications. In this manuscript, we, therefore, present a DG-based numerical scheme for our gas hydrate model, extended with the NCP-based semi-smooth Newton solver to handle the inequality constraints related with the phase transitions. In Section 2, we present the methane hydrate model based on Gupta et al. (2020). In Section 3, we outline our numerical algorithm based on the DG discretization scheme summarized in Appendix A. Finally, in Section 4, we present our numerical results. First, we validate the numerical scheme and its implementation by considering a 1D scenario of burial-driven recycling that was analyzed in Gupta et al. (2020), and second we simulate synthetic 2D scenarios of burial-driven recycling with different configurations of anomalous anisotropic high-permeability layers in the path of upward migrating gas. With these numerical examples, we demonstrate the numerical capabilities of our DG-scheme, and also highlight the necessity of modeling such layers with precision and care.

Mathematical Model
We consider the gas hydrate model developed by Gupta et al. (2020). The representative elementary volume (REV) for the model is shown in Figure 1, (Gupta et al., 2020).
This model is based on the theory of porous media and accounts for the following multiphysics processes: a) Dynamic evolution of the GHSZ due to changes in thermodynamic pressure-Temperature-Salinity (pTS) states, b) Migration of dissolved and gaseous methane through evolving GHSZ, c) Rate-based gas hydrate phase transitions, d) Exsolution-dissolution of methane in pore-water, and associated appearance and disappearance of free gas phase, e) Thermal effects, including the heat of hydrate phase change, f) Salinity changes, including feedbacks on methane solubility as well as hydrate stability, g) Changing sediment properties due to changes in pore-voids due to hydrate phase changes.
Moreover, we assume that salt is not allowed to precipitate as solid phase, sediment porosity is constant, that is, sediment compaction is neglected. A detailed model description including underlying assumptions can be found in Gupta et al. (2020). In the following, a summary of the main governing and constitutive equations is presented. PEIRAVIMINAEI ET AL. Let Ω be a bounded domain which contains the porous medium where the methane hydrate recycling (MHR) occurs. We assume that Ω is a subset of ℝ , ∈ {1, 2} and has smooth boundary (Lipschitz boundary). ∂Ω D and ∂Ω N denote the Dirichlet and Neumann parts of the boundary respectively and ∶= [0, ] denotes the time interval. Subscript β ∈ {h, g, w} denotes three pore-filling phases (gas(g), water(w) and hydrate (h)), subscript α ∈ {g, w} denotes two fluid phases (gas(g), water(w)), subscript s denotes the solid phase (sediment matrix), and superscript κ ∈ {M, H, c} denotes three components (Methane(M), Water(H) and salts(c)) in a porous medium. We introduce the following functions as model variables, Let U be the vector of model primary variables, which is some subset of the above introduced functions. We consider from the set of all model variables (Equation 1), as the primary variables.

Mass and Momentum Balance
Mass and momentum balance equations in a porous media for the components κ can be written as follows: where ρ α and μ α are the density and dynamic viscosity of the fluid phase α with velocity v α relative to primary sediment matrix. is diffusion flux of the component κ through fluid phase α. The time derivative terms accounts for the rate of change of the total mass of the component κ in each fluid phase ( ) . The second and the third give the amount of advection and diffusion of and is the volumetric source resulting from the hydrate phase change. Furthermore, we assume that there is no salt in the gas phase, i.e., = 0 .
We assume that hydrate phase is immobile, that is, v h = 0. Mass balance equation for the hydrate phase is given by Darcy's Law for the momentum balance of the fluid phase will be considered Ficks Law for the diffusive mass flux through the composite sediment matrix will be considered where k rα and are relative permeability and molecular diffusion coefficient of the component κ through fluid phase α respectively. τ is the sediment tortuosity. Since = 0 then = 0 . In addition, we have ∑ = 0.

Energy Balance
Energy balance equation is given by where ̇ℎ denotes heat of hydrate phase change and h α is the specific enthalpy of fluid phase α, u γ is the specific internal energy of the phase γ = g, w, h, s, and is the effective thermal conductivity,

Hydrate Phase Change Kinetics
We represent the rate of methane hydrate reaction, for a more detailed description see (Gupta et al., 2020). When methane hydrate is exposed to different effects such as depressurization or warming up, it decomposes to it's components. The rate of this reaction is modeled by the Kim-Bishnoi kinetic model (Kim et al., 1987). In this model, the rate of water and gaseous methane generated as a result of hydrate phase change are evaluated bẏ where p e is the equilibrium pressure for the methane hydrate, k r is the intrinsic reaction rate of hydrate phase change, and A rs is the reaction surface area available for hydrate phase change. is the molar mass of the component κ and N h denotes the hydration number. If ( − ) > 0 , the hydrate becomes unstable and If ( − ) < 0 , the hydrate becomes stable. Hydrate stability means that < 0 and according to mass balance Equation 3, if gaseous methane exists, it will be trapped between the water molecules and methane hydrate forms (Hydrate formation). Conversely instability means that > 0 and if methane hydrate exits, it releases gaseous methane and water molecules (Hydrate dissociation). ℎ denotes the hydrate formation (dissociation) rate. In addition, the following condition holdṡ+̇+̇ℎ = 0. where ℎ and the heat of hydrate phase change is given bẏ

Closure Relationships
The capillary pressure occurs across the gaseous and aqueous phase interface due to balancing of cohesive forces within the liquid and the adhesive forces between the liquid and soil matrix and relates the water and gas pressure.
Moreover, we have the following relation between phase saturations.
that is, total pore space is filled with hydrate, aqueous, and gaseous phases.

Vapor-Liquid Equilibrium (VLE)
Methane and water components are assumed to exist in vapor-liquid equilibrium and Henry's law and Raoult's law are valid. These equations relate mole fractions of each component existing in different phases.
Henry's law: where z M is the compressibility factor for methane gas, evaluated using Peng-Robinson equation of state, is the Henry's solubility coefficient for methane gas in water, and is the saturated vapor pressure for water in contact with methane gas. In addition to Equations 15 and 16, total mole fraction of components within each phase is bounded by 0 and 1, and it is 1 if the phase α is present, that is, The second inequality in Equation 17 holds when the corresponding phase disappears, that is, s α = 0, and because the VLE assumed to be valid, the mole fractions of the components remain undersaturated and are calculated from VLE Equations 15 and 16. We can write an equivalent complementarity conditions to the conditions in Equation 17 Now we consider the nonlinear complementary problem for the inequality constraints (Equation 18):

Hydraulic Properties
Capillary pressure is parameterized using Brooks-Corey model (Brooks & Corey, 1964) and is given by denotes the normalized aqueous phase saturation and S wr and Sgr are the irreducible aqueous and gaseous phase saturations respectively. Moreover, P 0 is the capillary entry pressure, λ is the material parameter related to sediment grain-size distribution, and m is the material parameter related to sphericity of hydrate growth.
Following the Brooks-Corey model, we parameterize the relative fluid phase permeabilities, Furthermore, we consider the effect of changing effective pore space due to hydrate phase change for the sediment permeability, where K 0 is the absolute permeability tensor of the primary sediment matrix. The exact functional relationships can be found in detail in (Gupta et al., 2020). H, h, c, e, ncp1, ncp2} be the set of indices of the corresponding Equations 3, 4, 8, 19, and 20. We remark that the number of equations is equal to the number of primary variables and thus after discretization we obtain a quadratic matrix. Moreover, let Ω ⊆ Ω and Ω ⊆ Ω be the corresponding Dirichlet and Neumann boundary conditions for U i .
Then, we have the following nonlinear problem, where U 0 , and are given functions.
Note that this problem is composed of a strongly coupled and highly nonlinear system of differential algebraic system of equations with four partial differential Equations 3, 8, one ordinary differential Equation 4, and two algebraic constraints (Equations 19-20).

Numerical Algorithm
Problem 1 is discretized in space using a DG method of order q defined on a quadrilateral mesh with ℎ elements and mesh size h. A fully Implicit Euler (IE) method is used to discretize the ODE system resulting from the spatial DG discretization. A brief description of the discretization scheme is given in Appendix A. The resulting nonlinear residual equations can be represented in compact form as follows, (Equation A22): where is the solution vector at time t n .
The nonlinear system (Equation 25) is linearized using a semi-smooth Newton solver, see (Wohlmuth, 2011) and references therein, which ensures that the rapidly switching phase states due to phase transitions remain consistent within each Newton step. In our scheme, we solve Equation 25 monolithically, that is, phase states will be determined along with solving the mass and energy balance equations. Within Newton loops, each phase state is determined by NCP Equations 19 and 20 which partitions the degrees of freedom (Dof) into active/inactive sets.
10.1029/2022JB025592 8 of 26 These active/inactive sets may change during the Newton loop, but the convergence of the Newton method guarantees the physically correct phase state of the system. The classical Newton method is applicable in the subsets of ℝ ℎ where the functional  is differentiable. Algebraic Equations 19 and 20 are semi-smooth and piecewise differentiable. We extend the values of the derivatives from nondifferentiable to differentiable regions using central difference method to calculate the Jacobian for our Newton scheme. The Jacobian can be numerically calculated using directional derivative: the ith column of the Jacobian can be obtained by where ϵ > 0 is a small positive number and ∈ ℝ ℎ , = 1, .., ℎ is the standard basis vector.
Time step sizes Δt n are adaptively adjusted according to heuristical rules based on the Newton performance, that is, Δt n increases by 20% if the Newton method converges in less than 5 iterations, decreases by 20% if the Newton method converges in more than 8 iterations, remains unchanged if the number of iterations of the Newton method is between 5 and 8, compared to Δt n−1 . If the Newton Solver does not converge within 10 iterations, then we redo the time step with a step length of Δ 2 . In the following, we outline our numerical algorithm 1 for solving Equation 25. 0 is the initial solution vector, j is the Newton iteration superscript, and  ( ) is the Jacobian matrix of the residual vector  at . Moreover, the algorithm contains the following numerical parameters: Δt max , Δt min → maximum and minimum time-steps allowed for the adaptive time-stepping; tol → maximum error accepted for the residual functional  ; j max → maximum number of Newton steps considered at each time-step; j 1 , j 2 → number of newton steps to adapt time-step sizes. Input: Δt max > Δt 0 > Δt min > 0, tol > 0, max, 1, 2 ∈ ℕ, 0 Output: 10.1029/2022JB025592 9 of 26 The numerical algorithm 1 is implemented using the software framework DUNE-PDElab (Bastian et al., 2010), version 2.7.0 (https://www.dune-project.org/modules/dune-pdelab/). To solve the linear system in line (1) of the algorithm, we use an in-built biconjugate gradient stabilized method (BiCGSTAB) (Blatt et al., 2016), as iterative solver for the linearized system. In all numerical examples presented in this manuscript, the parameters of the DG scheme in Appendix A and the numerical algorithm are chosen as max = 10, 1 = 4, 2 = 8, = 10 −6 , where σ is the penalty coefficient and Θ chooses the type of DG scheme (Θ = 1 for symmetric DG scheme, Θ = 0 for incomplete DG scheme, Θ = −1 for non-symmetric DG scheme).

Numerical Results
Our main motivation for the development of this new simulation framework based on the DG method arose from the need for accurate and robust handling of the multiphysics dynamics of the MHR problems in complex geological media, especially in relation with large local anisotropy and material heterogeneities. Our existing simulation environment (Gupta et al., 2020) is based on a finite volume based numerical scheme which offers many advantages like being fully locally mass-conservative, monotonic (i.e., no overshoots and undershoots, even with coarse mesh), conceptually simpler (in terms of implementation of active/in-active sets related with the semi-smooth Newton method), and computationally cheaper (due to low order and therefore fewer degrees of freedom). However, it has a major limitation when the subsurface properties show large local anisotropies and other complex material properties such as cross-cutting features like fractures and brecciated layers. With that in mind, we present here two numerical examples: (a) Example 1 considers a simplified 1D MHR scenario in a fully homogeneous medium with continuous burial at a constant rate. This example is used as a benchmark to validate the implementation of our DG scheme against our FV simulator, and (b) Example 2 simulates MHR scenarios in a more complex 2D setting where two different configurations of idealized anomalous anisotropic material layers are considered in the GHSZ. The goal of Example 2 is to demonstrate the capability of our simulator in handling such complex sediment structures, and to highlight the impacts on prediction accuracy that can arise from incomplete and/or inaccurate approximation of the properties of these complex sediment structures. Material properties are given in Tables 3 and 4. For our simulations, the domain of interest remains constant, despite the sedimentation, that is, at time t > 0 top of our domain of interest is not sea floor and it shifts downwards, red box in Figure 2, for example, at time t = 90 Kyr, sedimentation depth is 90 m (90,000 × 0.001), therefore, top boundary of the domain of interest lies at the depth of 90 m. In Figures 3, 4, and 6, we don't add sedimentation depth, because the goal is to see the evolution of the MHR process. Moreover, we consider hydrostatic pressure and local thermal equilibrium in both Examples inside the domain of interest.

Example 1. Validation scenario: MHR in a homogeneous domain
This scenario, developed and analyzed by Gupta et al. (2020), is based on the geological setting of a buried channel-levee (BCL) complex in the Danube paleo delta (Black Sea) that is believed to have deposited its levees between 320 and 75 kilo-annum before present (ka BP) (Zander et al., 2017). Here, we simulate how a continuous deposition of sediment layers over the past 300 ka could have affected the MHR though the gas hydrate stability zone (GHSZ). Hence, the initial setting is based on the paleo conditions existing at 300 ka BP, and the top of the computational domain is pinned at the corresponding paleo sea floor. We consider a 1D domain Ω = [ − 500, 0]. The problem schematic is shown in Figure 2.
Initial and boundary conditions are specified in Table 1. At t = 0, we assume hydrostatic pressure at sea floor p w = 15 MPa and a bottom water temperature of T = 4° C. We assumed that the initial pressure distribution within the computational domain follows a hydrostatic gradient, and the initial temperature distribution follows a steady-state geothermal gradient of 35° C/km. Moreover, water salinity is 3.5%. For detailed description of the initial and boundary conditions see (Gupta et al., 2020). For the prescribed paleo pTS conditions, the base of the GHSZ (bGHSZ) (i.e., depth at which p e = p g ) lies at 400 m below sea floor (mbsf). At t = 300 ka BP, we assume that there is no free gas anywhere in the domain, and methane hydrate is located in the interval [320, 400] mbsf, directly above the bGHSZ. Furthermore, we assume that the deposition of the sediment layer at z = 0 (i.e., paleo PEIRAVIMINAEI ET AL.

10.1029/2022JB025592
10 of 26 sea floor) occurs with constant sedimentation rate v s,z = 1 mm/year over a period of 300 ka (i.e., from 300 ka BP to present day).
As sedimentation occurs, more and more sediment layers accumulate above the paleo sea floor, leading to an increasing pressure and temperature at the paleo sea floor boundary. These changes cause the bGHSZ to shift  Gupta et al. (2020). left: The initial state of the system and the corresponding hydrate layer inside the gas hydrate stability zone (GHSZ), t = t 0 = 0 (i.e., 300 ka BP). right: The state of the system at t = t n > 0, indicating how the GHSZ shifts as a result of sedimentation over time. Note that the red box is the simulation domain for our 2D case. upwards and destabilize the overlying gas hydrate layer. As gas hydrates melt, methane is released, which in sufficiently high quantity can lead to a free-gas phase to form and accumulate at the base of the gas hydrate layer. Schmidt et al. (2022) have shown that gas migration through the GHSZ in this scenario is highly dynamic and occurs in cycles. The gas hydrate layer acts as a converging-diverging nozzle in the path of upward migrating free gas, Figure 5. To be consistent with the mathematical model, in Problem 1 we consider sedimentation time period of 300 ka from t = 0 to t = 300 Kyr where t = 0 corresponds to the 300 ka BP.
In Figure 3, to validate our scheme, we compare our results with the Example 1 from Gupta et al. (2020). Snapshots of methane hydrate, gas saturation, and salinity at t = 90, 150 Kyr for both FV scheme (Gupta et al., 2020), and our DG scheme are plotted.
Note. Initial salt mole fraction in aqueous phase is 0.0096, therefore, water salinity is 3.5%. We assume that initially there is no gas in our simulation domain, s g = 0, and dissolved methane = 0 . Moreover, hydrate layer lies between depth of 320 and 400 (m) with maximum saturation of 30%. Water mole fraction in gaseous phase is then determined by VLE equation. In Figure 4, snapshots of methane hydrate dissociation, gas migration, hydrate reformation is shown from 300 ka BP to present day. The process of gas migration through GHSZ and MHR can be summarized as follows: Free gas phase appears below the melting gas hydrate layer and flows upwards due to buoyancy. However, as the gas flows through the hydrate layer, the gas velocity continuously decreases because of decreasing permeability due to increasing hydrate saturation (converging part of the hydrate nozzle, see Figure 5). Once the gas phase passes the point with maximum hydrate saturation s h (throat of the hydrate nozzle, see Figure 5), the gas velocity starts to increase (diverging part of the hydrate nozzle, see Figure 5). As gas escapes the hydrate layer into the overlying GHSZ, reformation of methane hydrate occurs. The new hydrate layer continuously grows consuming the free gas provided by the dissociation of the previous methane hydrate layer, as shown in Figure 4.
While the increase of temperature by the geothermal gradient has a global effect on the equilibrium pressure, salinity has a local effect on the equilibrium pressure. This is due to the fact that heat diffuses much faster than salt (see Figure 6). Figure 7 shows linear to fourth order approximation of the solution of Example 1. The reference solution obtained on a fine mesh, h = 0.25 (m) with order q = 4 is plotted with the dashed line. To reduce the spurious oscillations in the gas saturation, we implemented a linear polynomial reconstruction (slope limiter) for the gas saturation based on the weighted mean derivatives of the solution in the neighboring elements (Frerichs & John, 2021). After convergence of Newton method at each time step, we restrict the slope of the approximated solution. If the absolute value of the slope of the approximated solution is bigger than the numerical central difference of the neighboring cell averages. then the approximated solution will be replaced by linear approximation based on neighboring cell averages, Figure 8a. In Figure 8a, snapshot of gas saturation is plotted for h = 0.5 and q = 1 at t = 100 Kyr with and without slope limiter. We chose t = 100 Kyr, because as it is shown in Figure 4, the maximum gas saturation occurs around t = 100 Kyr and the escaping of gas from hydrate layer leads to oscillatory sharp front. It shows that sharp gradients of gas saturation are avoided by implementing the slope limiter.
In Figure 9, the convergence behavior of the nonlinear solver for a mesh with h = 0.125(m) is shown. Time step sizes drop when gas phase appears. However, they mostly remain bigger than 10 years, while in (Gupta et al., 2020) Figure 4b, Gupta et al., compared the time step size of NCP and PVS approaches which showed that even for the coarser mesh, that is, h = 0.3125(m), the maximum time step size of 10 years was scarcely achieved.

Example 2. 2D scenario: Gas flow through GHSZ with heterogeneous material property
Here, we extend the above 1D scenario by introducing an anomalous material layer with high-permeability and large anisotropy within the paleo GHSZ. In the geological setting these layers represent pipe or chimney like

10.1029/2022JB025592
15 of 26 structure with bericciated sediment or anomalous sand lenz. We consider homogeneous sediment for background material to highlight the effect of the existence of such anomalies on the MHR process. The problem schematic is shown in Figure 10, and the initial and boundary conditions are similar to Example 1 and for completeness are given in Table 2. Two different configurations are considered for the permeability tensors of the anomalous material layer, as shown in Figure 10, where K 2 is permeability tensor of the anomalous layer that is rotation of background permeability tensor K 1 with θ degree and scaled by K F . Figure 9. Numerical results for Example 1, the evolution of the time-step size during the simulation, q = 1, h = 0.125(m). Time step sizes drop when gas phase appears, however, they mostly remain bigger than 10 years. Figure 10. Schematic setting of Example 2 showing the anomalous anisotropic layer. v g is the gas velocity, K 1 is the background permeability tensor and K 2 is the permeability tensor of the anomalous layer. In the geological setting these layers represent pipe or chimney like structure with bericciated sediment or anomalous sand lenz. Note that for both cases, we consider the scaling factor, K F , to be 100, that is, the anomalous anisotropic layer is 100 times more permeable than the background sediment.

10.1029/2022JB025592
16 of 26 K 0 is the absolute scalar permeability of the background sediment and K F is a scaling factor for the absolute permeability of the anomalous layer.
In the first configuration, the degree of rotation is 0°. This is the simplest form of anisotropy, and the only form that can be handled by our finite volume based numerical solver with linear two-point flux approximations. This form, however, ignores the strongly directional properties of such a layer (e.g., flow through fractures). The second configuration accounts for this by rotating the permeability tensor along the layer axis. Such form of anisotropy (with full tensor) cannot be handled with a linear two-point flux approximation in a finite volume scheme, and instead, requires advanced methods like multi-point flux approximations or non-linear two-point flux approximations, both of which are computationally more expensive and conceptually more complicated. However, in DG discretization, any form of material anisotropy can be handled easily without additional overheads. The direction of the gas velocity depends on the permeability tensor of the layers, as schematically shown in Figure 10. When gas reaches the high permeable layer in case of the rotated permeability tensor, it will flow dominantly along the layer, bypassing the regions above the layer, as shown in Figure 10b. Figure 11 shows snapshots of the MHR and gas migration processes for configuration with θ = 0°, and Figure 12 shows snapshots for the configuration with θ = 45°. Notice how in the former, more gas is transported to the region above the anomalous layer and a thicker hydrate layer with higher saturation develops, while in the latter, gas is completely diverted through the anomalous layer, and the gas migration through the GHSZ is fully localized within a focused flow channel that looks strikingly similar to the chimney-like fluid escape structures observed in seimic profiles, see (Crutchley et al., 2021;Waage et al., 2019). Notice also that the gas ascent toward the sea floor is much faster in the latter configuration compared to the former. What is especially interesting is that even though the anomalous layer has the same geometry and same heterogeneity in both configurations, the gas migration shows completely different behavior due to the nature of the anisotropy. These idealized scenarios clearly demonstrate how the approximation of the properties of the complex sediment structures can lead to remarkably large deviations in the system dynamics.
In fact, the development of focused gas flow in Figure 12c is a particularly interesting result with direct implications for real world scenarios. For example, seismic data of the gas hydrate system in New Zealand's southern Hikurangi subduction margin shows a network of normal faults lying within the GHSZ (Crutchley et al., 2021). Data shows a broad zone of both negative-and positive-polarity reflections (interpreted as sediment layer with coexisting free gas and gas hydrate). This zone lies directly beneath sub-vertical gas-flow conduits (possibly a combination of gas-charged normal fault and gas pipe/chimney), and extends up to the base of the regional GHSZ. In their analysis of this data, Crutchley et al. (2021) highlight the importance of considering the structural heterogeneity within GHSZ and in particular, the impact of normal faults on the gas migration through the GHSZ. Similarly, high resolution 3D seismic data from the Storfjordrenna gas hydrate pingos field in northwestern Barents Sea shows that the pingos lie on top of gas chimneys that are connected to inclined faults within the underlying free gas and hydrate?bearing sedimentary rocks (Waage et al., 2019), highlighting once again the relationship between gas hydrate dynamics and regional fault system. Our numerical scheme capture this dynamics quite well and results illustrate the role of structural heterogeneity on dynamics of gas hydrates and gas migration through GHSZ. Moreover, our results emphasize the dramatic deviations that can appear in simulated system behavior if the structural heterogeneities are not appropriately handled.

Conclusion
Natural gas hydrate systems are characterized by strongly coupled and highly dynamic multiphysics interactions that require sophisticated numerical schemes to capture the system behavior accurately and robustly. A particular challenge is related to the complex structure of the geological subsurface. Classically, problems like burial driven Boundary conditions x ∈ ∂Ω, t > 0 z = z sf p w = 15 MPa + ρ w gv s,z (t n + Δt) recycling are studied in homogeneous sediments, or sediments with a layered stratigraphy that follows the paleo and present sea floor topographies. Existence of strongly anisotropic anomalous layers with large contrasts in properties within the gas hydrate stability zone can lead to significant deviations in the system dynamics. An accurate prediction of these deviations is critical for estimating the present day geological carbon repositories, response of hydrate-bearing sediments to changing environmental and climate stressors, and their geomechanical stability in response to natural and anthropogenic activities. In this manuscript, we present a new numerical scheme based on the DG method for our methane hydrate model (Gupta et al., 2020). The motivation for the development of this new numerical scheme was to enhance the flexibility compared to FV approaches so that we can handle the structural complexities of the sediments more accurately, and therefore, be able to consider more realistic geological settings. The choice of the DG scheme was specifically inspired by the fact that it is locally mass conservative (like the finite volume method on which our earlier simulators are based), and can approximate the fluxes in anisotropic fields more generally without additional overheads (like larger stencils that are needed for extending finite volume schemes with methods like multi-point and nonlinear two-point flux approximations in order to capture material anisotropy). Here we show that (a) the semi-smooth Newton solver for handling gas-water phase transitions performs well with a DG based discretization, (b) the presented DG scheme is able to capture the multiphysics dynamics of the methane hydrate systems accurately, and (c) the presented DG scheme is able to accurately capture the gas migration and hydrate recycling processes through strongly anisotropic materials. We also demonstrate that layer properties influence sensitively the numerical simulation results and incomplete knowledge can result in very large prediction errors of the recycling process.   Table 4 Water and Methane Properties

Appendix A: Discontinuous Galerkin Method
In this section we present a discontinuous finite element method for Problem 1. In the following, all functions and parameters are assumed to be non-dimensional. The domain, Ω ⊂ ℝ , will be partitioned into quadrilateral elements Ω ∈ ℎ where ℎ is a mesh of the domain and = 1 … |ℎ| . Broken Sobolev spaces and Bochner spaces can therefore be written as F is called an interior interface if |F| ≠ 0 and there exist Ω − and Ω + in ℎ , see Figure A1, such that Ω − ∩ Ω + . n F is the unit normal vector to the interface F and the direction is arbitrary but fix. Let  be the set of all interior interfaces. Similarly, let  =  ∪  be the set of all the boundary faces including Dirichlet  and Neumann  boundary faces. Let ∈  and Ω − ∩ Ω + then ∀ ∈  1 (Ωℎ) , we introduce Moreover, the definitions for all boundary faces has to be adopted in a proper way. In the following, we consider flux continuity on the interfaces, that is, therefore, ∀ ∈  and ∀ ∈  1 (Ωℎ) , we have Let  =  ( ;  1 (Ω, ℎ) ) be the Bochner space for p w , T and  =  ( ;  1 (Ω, ℎ) ) be the Bochner space for , ℎ, , , , then ∀i ∈ C, we consider the following functionals to define variational formulation for Problem 1, Thus, variational formulation for Problem 1 can be written as follows Problem 2. (Variational Formula) Find ∈  2 ×  5 such that ∀ ∈  1 (Ω, ℎ) and ∀i ∈ C PEIRAVIMINAEI ET AL.

10.1029/2022JB025592
24 of 26 where ℎ = |ℎ| ( ) . The primary variables can be written in terms of the basis functions as follows: be the finite dimensional Bochner subspace for i ∈ {p w , T} and  ℎ =  ( ;  ℎ ) be the finite dimensional Bochner subspace for ∈ { , ℎ, , , } . Now interior and boundary penalty discontinuous Galerkin (IPDG) method for system (23) can be written as follows

A1. Implicit Euler Method
In this section, we apply Implicit Euler method for the system (A19). To do so, we consider a partition t 0 ≔ 0 < t 1 < … < t m ≔ t end of the time interval . By defining Δt n ≔ t n+1 − t n , We use finite difference approximation of the time derivative, where U n ≔ U (t n ). Let be the coefficient vector of the primary variables U w.r.t basis functions in  , that is, is a row vector with ℎ = ∑ ∈ ℎ elements.
Using finite difference approximation of the time derivative (A20), we introduce the residual functional  ( +1 , ) for every ∈  , = 1, … , ℎ and ∀i ∈ C as follows: The nonlinear residual equations can then be written in compact form as follows: Figure A1. Two neighboring cells with F as common interface.