Periodic states and their implications in gas hydrate systems

Highlights      •    Well-defined periodic states are embedded within the steady-state hydrate dynamics.     •    Periodic states lead to cyclic formation and dissociation of massive hydrate layers.     •    Periodic states are fully self-sustaining even in the absence of external triggers.     •    Spontaneous gas migration & pressure release occur in supposedly unperturbed systems.     •    Existence of periodic states implies an irreducible uncertainty in hydrate dynamics.  Abstract  Gas hydrates are one of the largest marine carbon reservoirs on Earth. The conventional understanding of hydrate dynamics assumes that the system, in the absence of external triggers, converges to a steady-state over geological time-scales, achieving fixed concentrations of gas hydrate and free gas phase. However, using a high-fidelity numerical model and consistently resolving phase states across multiple fluid-fluid and fluid-solid phase boundaries, we have identified well-defined periodic states embedded within hydrate system dynamics. These states lead to cyclic formation and dissolution of massive hydrate layers that is self-sustaining for the majority of natural marine settings. This previously unresolved characteristic could manifest as spontaneous gas migration and pressure release in, supposedly, unperturbed systems. Our findings show that the gas hydrate systems are not bound to have unique steady-state solutions. Instead, existence of periodic states introduces an irreducible, but, quantifiable uncertainty in gas hydrate dynamics which adds significant error bars to global gas hydrate inventory estimates.

section_content: Gas hydrates (GH) are ice-like crystalline solids formed by water and methane gas molecules under specific thermodynamic (i. e. pressure-temperate-salinity) conditions (Sloan and Koh, 2007). While uncertainty remains about their abundance in nature (Burwicz et al., 2011; Boswell and Collett, 2011; Milkov, 2004; Wallmann et al., 2012), they are commonly assumed to play an important role in Earth's carbon cycle with hydrate formation and dissociation being sources and sinks of methane (Wallmann et al., 2012). The scales and impacts of methane release from natural gas hydrate deposits have been investigated for various environments, especially those that are sensitive to climate change, such as the Arctic (James et al., 2016; Kretschmer et al., 2015), the continental margins and shelves (Johnson et al., 2019; Portilho-Ramos et al., 2018; Zhang et al., 2020), and the permafrost areas (Chuvilin et al., 2018; Frederick and Buffett, 2014). Accordingly, hydrates are often interpreted in terms of steady-states on geological time-scales, perturbed only by relatively rapid environmental changes (Wallmann et al., 2012). For instance, fluid and gas venting at the seafloor, seepage sites, and even submarine landslides are often associated with hydrate dissociation in response to an external trigger, such as warming bottom waters (Ketzer et al., 2020) or sea-level fluctuations (Cremiere et al., 2016; Wei et al., 2022). The observations of double or multiple bottom simulating reflectors (BSRs) are also commonly attributed to the external changes in sedimentation regimes (Zander et al., 2017) rather than internal gas hydrate system dynamics. In this view, a change in hydrate distribution is the consequence of bringing a supposedly stable hydrate system into an unstable (or out-of-equilibrium) thermodynamic state. This perspective is attractive, as it allows relating direct observations like methane seeps, pockmarks, or submarine landslides to changes in current or past environmental conditions. It also attributes a modulating role to hydrates in the global carbon cycle, e. g., in the form of positive feedback mechanisms such as hydrate melting under contemporaneous global warming conditions (Biastoch et al., 2011; Ruppel and Kessler, 2017) or due to depressurization upon post-glacial rebound (Wallmann et al., 2018). 

Here we explore a new complementary perspective by providing evidence for a natural, intrinsic periodicity in the gas hydrate dynamics manifested in the in-situ cyclic re-building and dissemination of massive gas hydrate and free gas volumes. Most strikingly, we show that natural hydrate systems exhibit this rich internal dynamics under constant environmental conditions, i. e. without any external environmental triggers. We resolve this hydrate system dynamics by using complex multi-physics interactions, where the formation and dissolution of dis-tinct gas hydrate layers is controlled by organic matter degradation, methanogenesis, continuous burial of phases, gas migration dynamics, and GH phase change kinetics within the gas hydrate stability zone (GHSZ) (Gupta et al., 2020; Wallmann et al., 2006 Wallmann et al.,, 2012;; Burwicz et al., 2011; Burwicz and Rupke, 2019). Conventional understanding, cultivated through broadly used simplified mathematical models, dictates that over geological time scales these processes converge to a steady state leading to stable gas hydrate and free gas concentrations. However, our new analyses of the well-established equations describing gas hydrate systems show that the long-term stability of the gas hydrate dynamical system is not quite as straightforward. By analyzing a large spectrum of geological parameters typical for marine settings, we have discovered that for the vast majority of hydrate-bearing geological settings, the steady-state solutions of GH dynamics exhibit stable internal periodicity (i. e. periodic growth and dissolution of massive gas hydrate layers coupled with spontaneous free gas migration through the GHSZ). Most importantly, the amplitude of change between the lowest and the highest GH saturations reached within one full cycle shows a remarkable variability across geological settings, and can reach significantly large values (e. g., 30 percent of the sediment pore space in the sampled parameter space in this study). In case of the standard mathematical models that do not resolve the cyclic states, the obtained solutions (i. e. fixed-steady-state gas hydrate and free gas saturations) will be bounded within the limits of the true cyclic solution. In that sense, the cyclic states can be seen as significant 'error bars' on the fixed-steady-state solutions. 

This discovery has profound implications as it challenges the conventional view that changes in hydrate distribution and the associated manifestations in observational data can only occur in response to external drivers, which implies that seafloor and sub-seafloor observations such as seafloor venting, cold seeps, pockmarks, and multiple BSRs do not need to be directly related to out-of-equilibrium subsurface conditions but could result from a self-sustaining internal cyclicity in the gas hydrate system dynamics. This discovery also points to systematic (i. e. irreducible) uncertainty embedded within high fidelity gas hydrate models, which has direct implications for the estimation of global carbon cycling, gas hydrate inventories, as well as the prediction of system responses to changing climate and environmental conditions on both short and long time scales.
section_name: Introduction

section_content: To analyze the behavior of the GH dynamical systems, we simulated a 1D burial-driven recycling problem based on a uniform geological setting with a broad spectrum of parameters that cover majority of the gas hydrate-bearing sediments globally. 

The model considers non-isothermal, compositional, multiphase flow of methane, porewater, and gas hydrates, and accounts for the following physical processes: 1) Advective flow of gas and water; 2) Capillary effects at the gas-water interface; 3) Burial of phases due to sediment deposition on the seafloor; 4) Hydrate phase changes across the gas-hydrate and dissolved CH 4 -hydrate phase boundaries; 5) Changes in permeability due to hydrate phase changes, 6) In-situ generation of methane through organic matter (OM) degradation via sulfate reduction, methanogenesis, and anaerobic oxidation of methane (AOM); 7) Pore-water salinity, transport of dissolved salts, fresh-water recycling, and its effect on hydrate phase stability, 8) Dissolution-exsolution of methane, and model degeneracy related to the localized appearance ↔ disappearance of the free-gas phase; and finally, 9) Thermal effects which arise due to the non-isothermal nature of the hydrate phase changes and the strong temperature dependence of the hydrate-gaswater phase equilibria. 

A detailed description of the governing equations, constitutive models, reaction network, model parameters and the numerical methodology is included in Appendix A. An in-depth technical description of the mathematical and numerical framework used in this study was previ- ously presented in (Gupta et al., 2020). A summary of all relevant model variables and parameters is provided in the Table 1 for an easy reference. 

The idealized 1D test setting is depicted in Fig. 1. The top-boundary is located at a depth = 0, and coincides with the latest seafloor at any instant of time. Continuous sedimentation buries all phases at a rate of. In general, a gas hydrate layer (e. g., pre-existing GH reservoir) remains in a stable state within the GHSZ. Under continuous sedimentation, this stable GH layer will be continuously buried with sediment grains below the base of the GHSZ, resulting in progressive hydrate dissociation and formation of a free gas phase. Buoyant gas will start to migrate upward to re-enter the GHSZ and crystallize in form of a new stable GH layer, ultimately enriching the pre-existing hydrate reservoir. This process, known as the gas hydrate re-cycling, has been previously described for natural gas hydrate systems as being characterized by high rates of sedimentation (Burwicz and Haeckel, 2020; Burwicz et al., 2017). The computational domain parameters, initial conditions, and boundary conditions are listed in Table 2. The total depth of the sediment column within our computational domain is chosen as = to cover the depth zones at which biogenic methane generation, phase transitions, and the build-up of free gas column take place in our model. At the start of the simulation, we assume that the sediment column is fully saturated in pore fluids of a given salinity, and no dissolved methane, free-gas, and hydrate layer are present anywhere in the column. The pore water pressure is assumed to be hydrostatic and the temperature distribution is in equilibrium with the regional thermal gradient Δ. At the top boundary, fixed bottom-water conditions are imposed, which translate to Dirichlet constraints for all primary variables (governed by PDEs), and at the bottom of the domain, Neumann fluxes are imposed. The growth of the sedimentary column starts with the deposition of additional sedimentary layers on top of the model ) marks the base of the theoretical gas hydrate stability zone (bGHSZ). Similarly, the intersection between the gas pressure ( ) and the phase equilibrium curve ( ). Three phase region where dissolved methane, gas hydrate, and free gas co-exist (given sufficient amount of dissolved CH 4 ) is marked as a dark blue field (i. e. III phase zone). In case there is no sufficient methane dissolved in pore fluids to overcome the free gas solubility limit, there are only two phases present in the pore space, i. e. dissolved methane fully saturating the pore fluid, and gas hydrate (light blue field and the part of the diagram marked with a crossed pattern). Consequently, in case pore fluids are not fully saturated in methane, gas hydrate and free gas phases are not stable (gray field in the plot marked as a I phase zone). 

domain according to the sedimentation rate defined for each modeling scenario. After deposition, the top of the sediments becomes the new seafloor and all sedimentary layers present in the model domain undergo steady-state compaction according to (Berner, 1980). In our previous work (Schmidt et al., 2022), we have isolated the modeling parameters that have the largest influence on the modeling outcomes. 

Overall, we ran all combinations of the parameters in the following ordered sets, resulting in a total of 297 simulation scenarios: 0. 11, 0. 12, 0. 13, 0. 14, 0. 15, 0. 16, 0. 17, 0. 18, 0. 185, 0. 19, 0. 195, 0. 2}. The chosen range of intrinsic permeability (10 -15 ≥ 0 [ 2 ] ≥ 10 -17 ) covers a wide range of values from coarse-grained sediments to hemipelagic clay, and is representative for gas hydrate-bearing settings (Wallmann et al., 2012). Since our goal is to analyze the gas migration and GH dynamics driven by burial processes, we primarily focus on the influence of the burial rate parameter on the steady-state system behavior. There are large uncertainties associated with bulk sedimentation rates over geological times. However, to highlight the internal cyclic dynamics, we use constant burial rates over the entire simulation period (i. e., all external forcings are eliminated to highlight the self-sustaining internal cyclic states). By analyzing multiple sedimentation regimes with burial rates in the range of 0. 04

we explore how the sedimentation rate lying within given uncertainty bounds affects the overall behavior of the GH dynamics.
section_name: Methodology

section_name: Results

section_content: First shown by Schmidt et al. (2022), the internal periodicity in gas hydrate system is a result of the so-called "hydrate nozzle effect" where the hydrate (solid) layer acts like a converging-diverging nozzle in the path of the upward migrating free gas, and the complex phase transitions between hydrate(solid) ↔ dissolved-gas(fluid) ↔ freegas(fluid) and evolving pore geometry due to hydrate(solid-to-fluid) phase changes lead to a "pulsating" gas flow. This nozzle effect is not a consequence of any new model equations. Rather, our model uses the already well-established set of equations, and this effect is captured through a mathematically consistent resolution of the phase-states, leading to accurate phase transitions especially within the II-phase and III-phase zones (shown in Fig. 1 ). 

In Fig. 2, we show the steady-state solutions for the reference scenario characterized with parameters 0 = 10 -16 m 2, 0 = 10 -18 mol∕ ( m 2 ⋅ Pa ⋅ s ), and = 0. 05 cm/a. It is evident that under continuous burial, even without any external climate and/or environmental perturbations, the GH reservoir can exhibit a very rich dynamics localized in the vicinity of the base of the gas hydrate stability zone (bGHSZ). We can see that the pressure-temperature-salinity (p-T-s) states show a very large amplitude of oscillation in the vicinity of the bGHSZ, but the strength of these oscillations is sharply reduced upon propagation through the overlying GH layers. Among these, the oscillation of the gas phase pressure is the most consequential because of its saw-toothed character, where the pressure builds over a long time but is released over a very short time (almost instantaneously on a geological timescale as shown in Fig. 2E ). This near-instantaneous pressure release could be sufficient to trigger mechanical instabilities on slopes or initiate fractures in deeper sediments. 

The associated bulk quantities of interest (QoI), namely, 1) location of the bGHSZ (Fig. 3A ) defined as the lowest depth at which =, where is the gas phase pressure and the equilibrium pressure at the free-gas ↔ GH phase boundary, 2) depth-integrated volume of hydrate per squared-meter of the seafloor ( ℎ, Fig. 3B ), and 3) peak hydrate saturation within the gas hydrate layer ( ℎ,, Fig. 3C ) are also plotted. These steady-state solutions and the bulk QoI show a clear periodicity over time, suggesting that under continuous burial, the GHSZ contracts and dilates and the hydrate layer builds-up and melts ad infinitum, with a change of ∼60 m in the thickness of the GHSZ and ∼15% in the total hydrate volume over each periodic cycle. Moreover, even though the time-period of the periodic cycles is relatively large (∼220 ka), the hydrate build-up phase and the corresponding GHSZcontraction appears to occur much faster (over a period of ∼30-50 ka), compared to the hydrate-melting phase and the corresponding GHSZdilation (over ∼170-150 ka). 

The evolution of the hydrate distribution along the sediment column is resolved over one time-period in Fig. 3F to show the hydrate build-up and melting phases more clearly. The hydrate layer obstructs upward gas migration by reducing the effective permeability (see Eqn. (A. 7)). Due to the convex shape of the hydrate distribution along the sediment column (see ℎ curves in Fig. 3F ), the gas velocity, which is proportional to the permeability (s. t., ∝ ∶= 0

( 1 - ℎ ) ), does not reduce uniformly across the whole hydrate layer. Rather, the section of the GH layer below the peak saturation decelerates the upwardmigrating gas, while the section above accelerates the gas in a manner analogous to a converging-diverging mechanical nozzle. If the upward migrating gas manages to seep past the throat of this hydrate-nozzle (i. e. peak saturation) this gas is flung across the hydrate layer to the overlying GHSZ, where the gas quickly converts back into hydrate, building a new layer above the pre-existing (or old) GH layer. The continuous
section_name: Periodic states

section_content: Initial and boundary conditions for the simulation scenarios. burial of the old layer below the GHSZ continues to supply the gas for the build-up of the new layer. Once the old layer is completely consumed, the new layer stops growing. Burial pushes this new layer towards the bGHSZ and the melting phase starts. The melting continues until enough gas has been built-up for it to escape past the peak of the GH layer, at which point a new cycle begins. The dynamics of hydrate nozzle is tightly controlled by the rate of gas supply at the base of GHSZ (proportional to the rate of hydrate dissociation) and the rate of upward gas migration (combination of upward buoyancy and deceleration from the hydrate nozzle). If the rate of dissociation is too low, the gas supply will be insufficient for its escape past the throat of the nozzle. Alternatively, if the rate of upward gas migration is too low (e. g., due to low intrinsic permeability 0 or high exponent or extremely high rate of hydrate formation), again the gas may not be able to escape past the throat of the nozzle. In both cases, the periodic states will not occur and the system dynamics will converge to a fixed steady state. Due to strong coupling of the hydrate ↔ dissolved-gas ↔ free-gas phase transitions with the p-T-s conditions, the cyclic building and melting of the hydrate also impacts the state of the GHSZ. The intersection of the phase boundaries leads to a distinct three-phase zone within the GHSZ where methane co-exists in all three phases. This zone is defined as the region in the sediment column bounded from below by the curve Γ ∶= - = 0 and from above by the curve Γ ∶= 4
section_name: Table 2

section_content: - 4 = 0. The evolution of the three-phase zone is resolved in Figs. 3D-E (yellow zone, where Γ is the point of intersection of the blue and red curves in part D and Γ is the intersection of the blue and red dotted curves in part E). The bGHSZ remains relatively stationary as long as the old hydrate layer is in place. When the old layer is consumed beyond a critical volume, the bGHSZ starts to shift upwards towards the base of the new GH layer. An additional point to note is that the hydrate layer enters the limit of the cyclic much faster whereas the gas saturation below the bGHSZ accumulates over multiple cycles before it converges with the true cyclic-steady-state. This is shown in Fig. 3F where the gas reservoir builds over 60 cycles. Given that the near-instantaneous pressure releases in each cycle would lead to mechanical instabilities sooner or later, it is unlikely that in nature the gas reservoir would ever build up to this theoretical steady-state gas volume. 

Finally, for this reference scenario, we also show the phase plots in the appendix in Fig. A. 9. The phase plots show the trajectories of the system variables relative to each other. These plots are an important tool to visualize whether the system converges to a fixed-point (true steady state) or a limit-cycle (i. e. periodic steady state), or exhibits signs of chaos with arbitrary trajectories. All variables of this reference scenario converge on a limit cycle instead of a fixed-point, showing categorically that for this particular combination of parameters, the gas hydrate dynamical system does not reach a true steady state, but exhibits a periodic steady state. Moreover, this periodic state appears to be a stable attractor, and at least within the simulated parameter space and time range, it does not exhibit signs of chaos. This does not mean that the GH system does not have chaotic states overall; only that this scenario does not. For models that do not resolve these cyclic states, we postulate that their estimations of the fixed-steady-states will lie within the limit cycles of the analogous periodic-steady-states. In that sense, we can see these periodic solutions to set error bars over the fixedsteady-state solutions.
section_name: Initial conditions Boundary conditions

section_content: So far, we have looked at the solutions of one particular (i. e. reference) scenario, and established the existence of periodic states that result from rich internal dynamics (as opposed to any external forcing of perturbations of environmental conditions). In Fig. 4, we plot the maximum and the minimum values of the QoI ℎ, that occur within one periodic cycle for all scenarios with 0 = 10 -16 m 2 and 0 = 10 -18 mol∕ ( m 2 ⋅ Pa ⋅ s ). Here, three distinct steady-state behaviors can be identified: I) No gas hydrate is formed in the domain, II) gas hydrate layer shows cyclic building, burial, and re-building, oscillating between two distinct peak values, and III) gas hydrate layer recycles continuously, reaching a fixed steady state value. The steady state behavior shows an abrupt transition from type-I to type-II to type-III back to type-I. This sudden qualitative or topological change in the behavior of the dynamical system for a small smooth change in a parameter value is called a bifurcation. 

On a process level, gas hydrate dynamics is a sum total of many competing processes. The parameter space explored here focuses on the competition between gas flow through upward migration (controlled by the permeability 0 ) and burial (controlled by sediment velocity ), and flow modulation through the converging-diverging hydrate nozzle (controlled by the reaction rate 0 ). If burial velocity is too low, OM degradation remains confined to the upper sediment layers and methanogenesis does not occur. If burial velocity is too high, the dissolved methane gas is rapidly buried below the base of the GHSZ and is unable to accumulate to form a free gas phase. In both these cases, gas hydrate is not formed and the system dynamics exhibits a type-I steady state. If the burial velocity is large enough to transport OM deeper into the sediment to promote methanogenesis, but small enough to allow the accumulation of methane into a free gas phase, gas hydrate layer will build up and undergo recycling in the vicinity of the base of the GHSZ, and the system dynamics will exhibit either a type-II or a type-III steady state, depending on how the hydrate nozzle modulates the flow of the upward migrating free gas. 

In Fig. 5, the steady states of all 0 - 0 scenarios are plotted versus burial velocity. Results highlight that the bifurcation manifold of a gas hydrate system has a highly complex shape in a high dimensional parameter space. A comparison of the steady-states suggests that high permeability and high hydrate kinetic rates promote the existence of periodic steady-states. Moreover, the shape of the bifurcation manifold in the 0 - plane corresponding to 0 = 10 -18 (Fig. 5B, E ) is nonmonotonic (i. e., the periodic-states envelope expands from 0 = 10 -17 to 0 = 10 -16 but contracts from 0 = 10 -16 to 0 = 10 -15 ), suggesting that in the hydrate-nozzle dynamics, the parameters 0 and 0 are not completely independent, and there likely exists an optimal pair with the largest range of cyclic states. any periodic states, so time-period is not defined. Gas hydrate dynamics is a sum total of many competing processes. The parameter space explored here focuses on the competition between gas flow through upward migration (controlled by permeability 0 ) and burial (controlled by burial rate ) and flow modulation through the converging-diverging hydrate nozzle (controlled by reaction rate 0 ). 

Finally, the time periods of all cyclic solutions are plotted in Fig. 5G-I. Results suggest that a combination of high permeability, low hydrate kinetic rate, and low burial velocity leads to shorter cycles. 

The existence of bifurcations in gas hydrate dynamics is an important discovery for many practical reasons. For example, in the parameter space sampled within this study, periodic states do not occur everywhere. Rather, they occur only over a range of burial velocities. The amplitudes of the periodic states and their time periods also show a huge variability across the 0 and 0 dimensions. Therefore, gas hydrate dynamics for parameters located within the bifurcation manifold will have a significant systematic (i. e., irreducible) uncertainty. This uncertainty is irreducible in the sense that it is a fundamental mathematical property of the system, unlike the other forms of parameter uncertainties (e. g. due to measurement errors or incomplete sampling, etc. ) that can, at least in theory, be eliminated by improving measurements, sufficient data, etc.
section_name: Bifurcation

section_content: The GHSZ is a function of the p-T-s states, and in general, the gas hydrate dynamics is strongly influenced by the prevailing hydrostatic pressure and thermal gradient. To get an impres-sion of the abundance of the periodic states in nature, we further tested the burial-driven GH dynamics over a range of water-column depths ( = {2195, 1097. 5, 548. 75} m ) and regional thermal gradients ( Δ = {45, 34. 22, 24. 8} C/km

). The GHSZs for the corresponding scenarios are shown in Fig. 6. Note that = 2195 m and Δ = 45 C/km were used to define the reference test-setting in Table 2. To test the sensitivity of the GH dynamics to pressure conditions, we chose additional values of ∕2 and ∕4. The Δ values were chosen such that the diagonal elements in Fig. 6 have the same depth of the GHSZ. An additional difference between the reference scenarios (Table 2 ) and these scenarios is in the bottom water temperature, chosen here as | =0 = 0 C. This was done to ensure that the shallowest GHSZ extends below the seafloor. Fig. 7 shows the bifurcation plots for each -Δ scenario with respect to the evolution of the hydrate saturation amplitudes over the burial rates. Corresponding bifurcation plots for specific volume of hydrates, location of the base of the GHSZ, and time-period of the periodic states are shown in the appendix. 

A few important observations can be made:

i There is no apparent correlation between the periodic dynamics of GH and the depth of the GHSZ. The GH dynamics is significantly different for the same GHSZ depths (parts A,E,I in Figs. 7, A. 10, A. 11, A. 12). ii A combination of shallow water-column and low thermal gradient lead to the appearance of chaotic states, especially at lower burial rates (parts A,B,E in Figs. 7, A. 10, A. 11, A. 12). The chaotic states are inferred from the phase plots, shown in Figs. 7, A. 10, A. 11, A. 12, where the time evolution of the variables relative to one-another shows arbitrary trajectories. This is in contrast with the stable periodic states, which converge to a fixed closed-loop called the limit cycle. The chaotic states in shallow sediments (part A) show an average time-period of gas release in the order of 10-to-1 ka, which is several orders of magnitude (∼2-3) smaller than the analogous periodic states. iii In shallow sediments, instabilities related to the onset of periodic states can lead to a 'choking' of the whole sediment column with hydrates (part D in Figs. 7, A. 10, A. 11, A. 12). iv Shallow water-column depth and high geothermal gradient lead to lower hydrate volumes ( 2 ), it appears that the periodic states are more abundant at lower bottom water temperature, although a more systematic analysis of the temperature effects is necessary.
section_name: Effects of hydrostatic pressure and thermal gradient

section_content: The GH-dynamics is controlled by the hydrate nozzle effect where the change in hydrate permeability due to hydrate phase transitions modulates the gas flow, leading to the periodic states. It is well known that the permeability evolution with respect to the hydrate phase change strongly depends on whether the hydrate coats the sediment grains or fills the pores as it forms (Kleinberg et al., 2003). Permeability decreases more rapidly with hydrate formation for the pore-filling than for the grain-coating growth habit. In our model, this is captured by the parameter (Eqn. (A. 7)), as shown in Fig. 8A. The occurrence of periodic states cannot be attributed to some specific parameterization of the - ℎ relationship or its growth habit. Instead, it is the consequence of the 'process' of flow modulation by the hydrate layer. The growth habit, however, does affect the properties of the emergent periodic states, as shown in Fig. 8B-E. In general, grain-coating hydrates lead to higher GH volumes and saturations, while pore-filling hydrates show much higher variability and amplitudes of GH volume, saturation, and bGHSZ, especially at higher burial rates.
section_name: Effects of hydrate growth habit

section_content: Bifurcations and periodic solutions are known mathematical properties of non-linear dynamical systems, and are generally studied under the broader umbrella of the chaos theory (Awrejcewicz, 2012; Kaneko and Tsuda, 2011; Feldman, 2019). Bifurcations and periodic steadystate solutions are well-known from other applications, e. g. Navier-Stokes flow (e. g., Kato, 1997), climate dynamics (e. g., Ghil and Tavantzis, 1983; Korobeinikov and Mcnabb, 2001), chemical kinetics (e. g., Nielsen et al., 1991), geomorphodynamics (e. g., Goehring, 2013), etc. It has, however, never been observed, or at least reported, in gas hydrate models before due to certain assumptions and simplifications made in the previous numerical models that predict averaged solutions (bounded within the limit cycle, as shown in the phase plots in the manuscript). Fig. 7. Bifurcation plots showing the evolution of the amplitudes of hydrate saturation over the burial rates for the scenarios that test the influence of the pressure and thermal gradient on the GH-dynamics and abundance of the periodic states. The orange shaded regions in parts A, B, and E highlight the chaotic states and the blue shaded region in part D highlights the solutions where the instabilities at the onset of periodic states led to a 'choking' of the sediment column with GH. Also shown are selected phase-plots that show the relative trajectories of the variables in chaotic, choked, as well as stable periodic states. 

Commonly used multi-phase gas hydrate models are based on the same mass balance and energy equations shown here, coupled with similar reaction network including kinetically controlled in-situ POC degradation, methanogenesis, anaerobic methane oxidation (AOM), and sulfate reduction processes. Therefore, in this manuscript, we do not show a 'new modelling method' that can capture these oscillations. Rather, we show that to capture the internal gas hydrate system cyclicity, it is necessary to, 1) consistently resolve the phase transitions occurring simultaneously across multiple fluid-fluid and fluid-solid phase boundaries, and 2) fully resolve the fluid-sediment interactions (i. e. parametric function describing permeability evolution with hydrate saturation, necessary to simulate the nozzle effect). When the used formulation is mathematically or numerically smoothed or simplified, the numerical solution will not show the cyclic states but will lie within the bounded limit cycles shown in Fig. A. 9. 

The internal periodicity is shown to be a self-sustaining process regardless of external forcing factors, such as anthropogenic warming, sea-level fluctuation, or large scale planetary cycles e. g. the Milankovitch cycle. However, the hydrate system response to overlapping short-and long-time scale forcing is yet to be analyzed. For instance, there is an apparent time lag in development of full system steady-state cyclicity observed in both gas hydrate and free gas phases. The gas hydrate phase reaches a full steady-state cyclicity several cycles before a steady-state free gas cyclicity. It is explained by the fact that free gas phase needs sufficient time to build up enough gas volume that allows upward migration. This is directly linked to the fact that the free Fig. 8. A) Permeability -hydrate saturation ( - ℎ ) relationship depending on the hydrate growth habit, where the parameter '' controls its type; B) Bifurcation plot showing the evolution of the specific hydrate volume vs. the burial rates for each growth habit; C) Time periods of the cyclic states vs. the burial rates for each growth habit; D) Bifurcation plot showing the evolution of the hydrate saturation vs. the burial rates for each growth habit; and, E) Bifurcation plot showing the evolution of the base of the GHSZs (bGHSZ) vs. the burial rates for each growth habit. gas saturation is limited by sediment permeability within the free gas zone, whereas gas hydrate peak is limited by kinetic reaction rate of phase transition. As a consequence, full steady-state system cyclicity of both gas hydrate and free gas phases can be reached faster in settings characterized by: a) high organic matter content available for efficient biodegradation, b) active fluid flow or high-permeability pathways for gas transport from below, b) additional methane sources, for example, gas reservoirs. However, it has been shown that a steady-state gas hydrate system cyclicity can fully develop in geological settings without constant gas supply (Schmidt et al., 2022) as a result of in-situ organic matter decomposition but could be fueled by sufficient amounts of free gas present in the sediment pore space as a stationary gas pocket. We attribute this observation to the fact that gas hydrate system cyclicity by itself is a consequence of the GH-nozzle formation. However, the nozzle dynamics is modulated by the supply and connectivity of free-gas from the reservoir at its base. To this effect, biogeochemical reactions provide strong feedbacks to the GH-dynamics, but are not 'necessary' for the appearance of the periodic states. Shallow hydrate systems (i. e. defined by a relatively thin GHSZ) might manifest hydrate cyclicity in features like seafloor venting of fluid and gas, slope instability, formation of pipe structures within the GHSZ, seepage sites, or moving BSRs. In contrast, deep hydrate systems might not exhibit observable changes at the seafloor.
section_name: Discussion

section_content: The apparent existence of bifurcation manifolds and in particular, the periodic states, has profound implications as it sets hard limits on the predictability of present-day gas hydrates through steady-state analysis. What this means is, that if the parameters for a particular geological scenario lead to a cyclic (or periodic) steady-state, we cannot say with confidence which part of the cycle we are on at the present moment. However, what we can predict with some confidence is the maximum and minimum amounts of GH that can occur for the given set of parameters. This means that there is a systematic irreducible source of uncertainty embedded within the GH dynamics. Given the large variation in the GH saturation and volume occurring over one time-period of the cyclic solution (e. g. ℎ ranging from 30 vol. % up to 80 vol. % and S. Gupta, E. Burwicz-Galerne, C. Schmidt et al. ℎ from 10 0 to over 10 2 m 3 /m 2 across scenarios shown in this study), current estimates of the present-day global GH inventory may contain a large margin of error. 

The periodicity in the steady-state solutions of p-T-s states also means that in typical geological settings, the GH systems are highly dynamic even without any external climate and/or environmental perturbations. Especially in relatively shallow gas hydrate systems, the landslides, slope failures, pockmarks, pipes, and chimneys etc. observed may have occurred spontaneously without any external triggers like sea level, bottom-water, sediment loading, or salinity fluctuations. 

The cyclic rebuilding of GH layers occurs in phases through the hydrate nozzle mechanism, where two GH layers can coexist, one at the base of the current/latest configuration of the GHSZ and one below (corresponding to some past configuration of GHSZ). The dynamics of the cyclic states can, therefore, also explain some of the observed double (and even multiple) BSRs. 

The existence of periodic states has particularly serious implications for the prediction of future climate impacts due to modern-day climate dynamics. It is likely that the anthropogenic climate perturbations (i. e., very fast changes in climate conditions) may push the formerly stable steady state of the GH reservoirs (estimated based on paleo-climate conditions) to new periodic states with large p-T-s fluctuations. Since the pressure release at the end of the cycle is near instant on the geological time-scales and rather significant (e. g. up to several bars), it can trigger mechanical sediment failures, including the formation of pipes, chimneys, or pockmarks (Cartwright et al., 2021). Thus, if formerly stable states are pushed into periodic states, the risks of mechanical failures and uncontrolled gas release will increase. Moreover, the anthropogenic climate fluctuations may also push the formerly periodic states (with longer time periods) towards chaotic states (with several orders of magnitude smaller time-periods), making long term predictions nearly impossible. 

The discovery of the periodic states in GH dynamics has opened up multiple fundamental questions related to past, present, and future evolution of natural GH systems. Some of the properties of these states, like self-sustaining dynamics in the absence of external triggers, sudden pressure release at regular time-intervals, spontaneous escape of trapped gases from GH layer, and a continuously out-of-equilibrium GHSZ, are likely to have considerable implications for coastal dynamics, evolution of continental margins, and seafloor/sub-seafloor infrastructure. We explored these properties along some limited sections of the parameter-space, but from the shape of the emerging bifurcation manifold, it is evident that the GH dynamical system is high-dimensional and much more complex. We believe that a full interpretation of the implications of these properties and their role within the solid Earth systems will require more systematic mathematical analyses besides scenario analyses and numerical simulations. in the domain vs, and F) QoIs ℎ, vs ℎ. The phase plots show the trajectories of the system variables relative to each other. These plots show whether the system converges to a fixed-point (true steady state) or a closed-loop trajectory called a limit-cycle (i. e. periodic state), or exhibits signs of chaos with arbitrary trajectories. All solutions for this scenario converge to a limit cycle instead of a fixed-point, showing that for this particular combination of parameters, the gas hydrate dynamical system does not and will not reach a true steady-state. Instead, it exhibits a periodic steady-state which would contain a simplified steady-state solution within its envelope. flow through this reactive media, an effective porosity is defined, s. t.,

The sediment phase constitutes the primary matrix which is assumed to be rigid (i. e. undeformable), whereas, the hydrate and sediment phases together constitute the composite matrix. The void-spaces embedded within the primary matrix define the porosity of the medium, which remains constant over time (due to the rigidity assumption). The hydrates, however, are reactive solids that can undergo massive volume changes. Therefore, the effective porosity of the composite matrix evolves over time, depending on the progression of the hydrate phase transitions. The evolution of the effective porosity is strongly coupled with the reaction kinetics of hydrate phase change, and provides highly nonlinear and bidirectional feedbacks to fluid flow fields. The model also considers compositional flow, where gas, water, and sediment phases are composed of multiple reactive species that can interact and undergo phase transitions across the fluid-fluid and fluid-solid interfaces. However, the changes in volume of the sediment phase are assumed to be negligible compared to those of the hydrate phase, and therefore, the impacts of sediment phase transition on porosity evolution are ignored in this model.
section_name: Global implications

section_content: Broadly, the model accounts for the following physical processes:

-Advective flow of gas and water; -Capillary effects at the gas-water interface; -Burial of phases due to sediment deposition on the seafloor along the continental margins; -Hydrate phase changes due to precipitation ↔ dissolution of hydrates in equilibrium with dissolved methane, and formation ↔ melting of hydrate in equilibrium with methane in the free-gas phase;

-Changes in hydraulic properties (i. e., permeability, capillary entry pressure, specific surface area) due to evolving porosity as a result of hydrate phase changes, -In-situ generation of methane through organic matter (OM) degradation via sulfate reduction, methanogenesis, and anaerobic oxidation of methane (AOM); -pore-water salinity, transport of dissolved salts, fresh-water recycling, and its effect on hydrate phase stability, -Dissolution-exsolution of methane, and model degeneracy related to the localized appearance ↔ disappearance of the free-gas phase; and, -Thermal effects which arise due to the non-isothermal nature of the hydrate phase changes and the strong temperature dependence of the hydrate-gas-water phase equilibria. 

The main governing equations, derived from the mass, momentum, and energy conservation principles, are briefly outlined as: where, denotes the phase pressures related through a pressure jump across the phase interface, also called the capillary pressure s. t. - ∶=, denotes the gravity, Θ the phase-wise flow-mobility, and the absolute permeability of the composite matrix, s. t.,

with 0 as the absolute permeability of the primary matrix. Please note that in this study, we do not consider the effects of fracture permeability on gas hydrate growth which would be an important consideration in clay-rich sedimentary facies. Furthermore, the variable denotes the burial velocity related to the rate of sediment deposition on the seafloor, s. t., the total velocity of any phase undergoing burial is given as,

where, and are Darcy velocities and ℎ = 0. 

The terms denote non-linear reactive-sources for each phase corresponding to following phase transitions: gas (methane) dissolutionexsolution ( = ), hydrate formation-dissociation ( = ), and hydrate precipitation-dissolution ( = ). 

The kinetic rate of the hydrate precipitation-dissolution transition is,

where, ℎ is the hydrate solubility, and is the rate of precipitation-dissolution phase change, s. t.,

The related source terms are,

The kinetic rate of the hydrate formation-dissociation transitions is,

where,

where, is the hydrate equilibrium pressure, 0 is the surface area of the hydrate-free sediment, and is the intrinsic reaction rate, with where, 4 + ℎ 2 = ℎ Note that the source term for gas dissolution-exsolution does not appear explicitly in the governing equations, because it cancels out in the total mass balance for methane in dissolved and free-gas forms. The generalized compositional system considers number of species, partitioned across the water and sediment phases. The ordered set of all species is denoted with  with cardinality () =. The variable denotes the concentration of the ℎ species in mmol per litre pore-water volume, where ∈ , and Φ transforms the concentration of the ℎ species to mmol per litre total volume, s. t.,

where, [] × is a mobility matrix s. t.,  = 1 if = and corresponds to an w-species, otherwise  = 0.  is the transport operator of the ℎ species describing the convective and diffusive mass fluxes, s. t.,

where, and are the molecular diffusion coefficients in water and sediment phases, respectively. Furthermore, the terms ∑ =1 and ∑ =1 are the chemical sources where and denote the stoichiometric coefficients for the ℎ species in the ℎ equilibrium and kinetic reaction, respectively, and and denote the rates of the ℎ equilibrium and kinetic reaction, respectively. The equilibrium reaction rates are unknown and cannot be estimated apriori. We find a suitable matrix In this manuscript, the compositional system is based on the OM degradation model described by Wallmann et al. (2006). The reaction network, consisting of the sulfate reduction, methanogenesis, and AOM reactions are summarized in Table A. 3. The compositional system also includes chlorinity, which is used as a proxy for pore-water salinity. 

The compositional system excludes the free-gas phase because freegas is composed of only one component, methane. Under an assumption of vapor-liquid equilibrium, if the concentration of dissolved methane remains below its solubility limit 4 in pore-water, free-gas phase cannot exist. Methane released through methanogenesis and melting/dissolving gas hydrates dissolves into pore-water until solubility limit is reached, beyond which all excess methane is spontaneously exsolved into a newly appearing free-gas phase. Conversely, perturbation of local pressure-temperature-salinity (p-T-s) state may raise the solubility limit s. t. the methane in the free-gas phase dissolves back into pore-water, leading to locally disappearing free-gas phase. Mathematically, the transition of methane across dissolved and free-gas states can be described by the following set of inequalities, Together, these lead to the Kharush-Kuhn-Tucker (KKT) type algebraic constraint given in eqn. (A. 4). This constraint converts the model into a constrained optimization problem. Moreover, similar to the equilibrium reaction rates in the compositional system, the rate of mass-exchange across the gas-water interface due to methane dissolution-exsolution cannot be known apriori. Therefore, we eliminate these terms by summing up the respective phase-wise and component-wise mass balance equations, leading to the conservation of total methane given in eqn. (A. 3). The factor 4 appearing in eqn. (A. 3) is the molar mass of methane, used to make the units of density and concentration consistent. 

Finally, in the energy conservation eqn. (A. 6), the variable denotes the homogenized temperature, the phase-wise specific heat capacity at constant volume, the specific heat capacity of fluid-phases, ℎ the phase-wise thermal conductivity, and the heat sources associated with each fluid-fluid and fluid-solid phase-transition. Eqns. (A. 1)-(A. 6) together form a strongly coupled and highly nonlinear system of partial-differential-algebraic equations (PDAE), to be solved for the primary variables: ]
section_name: A.2. Governing equations

section_content: The numerical scheme is based on a fully upwinded cell-centered finite volumes method for spatial discretization and an implicit Euler method for temporal discretization. The scheme is implemented in DUNE-PDELab (version 2. 8) (Bastian et al., 2010) based on C++. For the linearization of the system of governing PDEs, we have implemented a semi-smooth Newton solver which can handle the gas-water phase transitions and appearing and disappearing free-gas phase in a mathematically consistent manner. A highly optimized SuperLU (Demmel et al., 1999) linear solver is used to perform 1D calculations in sequential mode. In general, the numerical implementation is capable of solving in 1D, 2D and 3D. The 2D and 3D calculations can be performed OpenMPI parallel mode using a built-in Algebraic Multi-Grid solver. 

The computations for this study were performed on the highperformance computing cluster at Kiel University (CAU). Further details of our numerical scheme can be found in (Gupta et al., 2020).
section_name: A.3. Numerical solution

section_name: Acknowledgements

section_content: We are grateful to Laurence Coogan, Tim Minshull and two anonymous reviewers for providing us valuable comments that greatly improved our manuscript. SG acknowledges additional support from the SMART Project funded through the Helmholtz European Partnering Initiative (Project ID Number PIE-0004 ) and the MARCAN Project funded through the European Research Council (Grant No 677898 ) under European Union's Horizon 2020 research program.
section_name: Acknowledgements

section_name: Funding

section_content: This research is funded by the MSCA Postdoctoral ERA Fellowships 2021 action, under the Horizon Europe program, project ' WarmArctic ', number 101090338 and by the Cluster of Excellence ' The Ocean Floor -Earth's Uncharted Interface ' ( EXC 2077 ) funded by Deutsche Forschungsgemeinschaft (DFG) -Project number 390741603 hosted by the Research Faculty MARUM of the University of Bremen, Germany.
section_name: Funding

section_name: Data availability

section_content: The source code for the model and test scenarios is publicly archived and can be accessed as https://doi. org /10. 5281 /zenodo. 7920830.
section_name: Data availability

section_name: Data and materials availability

section_content: The numerical model is implemented within version 2. 8 of C++ based DUNE-PDELab framework (Bastian et al., 2010; Sander, 2020), and uses the in-built matrix assembler, linearization algorithm (Newton method with numerical Jacobian), and linear solver (parallel Algebraic Multi-Grid (AMG) solver with stabilized bi-CG preconditioner). The DUNE libraries used in this study are preserved at https://gitlab. dune -project. org /pdelab /dune -pdelab and are developed openly at https:// www. dune -project. org/. 

The source code for the model and test scenarios presented in this manuscript is publicly archived and can be accessed as https://doi. org /10. 5281 /zenodo. 7920830. Detailed step-by-step instructions are included in this repository (README. md) on how to install and run the code. Also included are numerical solutions for the reference scenario, along with matlab/octave based postprocessing routines to visualize these solutions.
section_name: Data and materials availability

section_name: CRediT authorship contribution statement

section_content: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
section_name: Declaration of competing interest

section_content: Here, we present a detailed description of the governing equations, reaction network, constitutive models, and model parameters used to simulate the burial-driven gas hydrate dynamics presented in this manuscript.
section_name: Appendix A. Mathematical model

section_content: The continuum-scale mathematical description of the conservation laws is based on the following homogenized variables defined over an REV (representative elementary volume (Helmig, 1997) Furthermore, the gas and water phases are mobile fluids, while the hydrate phase is an immobile solid that is chemically active, undergoing volume changes within the pores. Therefore, to describe the fluid Appendix B. Supplementary material An animation of the numerical solutions for the reference scenario (i. e., the scenario with permeability 0 = 10 -16 m 2, kinetic rate 0 = 10 -18 mol∕( 2 ⋅ Pa ⋅ ) and burial rate = 0. 05 cm/a) is included as supplementary information. The public git repository for the code includes the raw numerical solution files for this numerical simulation. 

Supplementary material related to this article can be found online at https://doi. org /10. 1016 /j. epsl. 2023. 118445.
section_name: A.1. Preliminaries