The Response of Gas Hydrates to Tectonic Uplift

Pressure reduction following uplift may lead to dissociation of gas hydrates. The dynamics of hydrate dissociation in such settings, however, are poorly understood. We used TOUGH+HYDRATE to investigate the response of gas hydrates to an uplift of 0.009 myr-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document} over the last 8 kyrs, the approximate end of the postglacial sea-level rise. Geological parameters for the simulations are based on hydrate deposits from the Nankai Trough subduction zone. Our results suggest stabilisation from endothermic cooling, elevated pore pressure, and pore water freshening significantly slows hydrate dissociation such that the hydrate remains in place at its pre-uplift level. A shallower hydrate layer forms from upward-migrating gas when assuming moderate to high permeability (10-15\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-15}$$\end{document} and 10-13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-13}$$\end{document} m2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2}$$\end{document}), while gas remains trapped for low permeability (10-17\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-17}$$\end{document} m2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2}$$\end{document}). In the latter case, we predict elevated pore pressure with potential implications for seafloor stability. Our findings suggest that following uplift, hydrates may exist outside the predicted regional gas hydrate stability field for thousands of years.


Introduction
Gas hydrates are naturally occurring ice-like solids, rich in carbon, e.g. methane gas (Sloan and Koh 2007). Methane hydrates sourced from biogenic methane gas are the commonest form of gas hydrates and are commonly found in high-pressure and low-temperature regions such as the permafrost and oceanic sediment pore spaces of continental margins (Kvenvolden and Lorenson 2001;Hester and Brewer 2009). According to Milkov (2004) and Boswell and Collett (2011), the estimated amount of carbon found in methane hydrates exceeds that of conventional fossil fuel accumulations, potentially making gas hydrates an energy source. However, a change in the thermodynamic conditions for gas hydrate stability could lead to the release of large amounts of methane, thereby affecting seafloor stability, global carbon cycle, and climate (Dickens 2001;Nisbet and Chappellaz 2009;Dessandier et al. 2020).
The base of the gas hydrate stability zone (BHSZ) is often marked by bottom simulating reflections (BSRs) in seismic profiles, usually caused by the presence of free gas (Singh et al. 1993;Rajput and Thakur 2016). The BHSZ constitutes the depth at which the geothermal gradient and pressure profile cross the gas hydrate phase boundary and thus should only occur at a single depth beneath the seafloor. However, double or multiple BSR(s) have been interpreted at several locations, e.g. Nankai Trough (Foucher et al. 2002). Furthermore, it is generally thought that double or multiple BSR(s) form when gas remain trapped at the paleo-BHSZ because of the low relative permeability of fine-grained sediments (Judd and Hovland 2009). They could also be associated with thermogenic gas, i.e. gas hydrates containing methane gas and significant quantities of higher hydrocarbons, e.g. propane (Geletti and Busetti 2011) and diagenetic positive-polarity BSRs (Berndt et al. 2004) at some locations. In addition, gas hydrates related to double or multiple BSR(s) have also been linked to geological processes associated with pressure-temperature changes, e.g. glacial cycle (Bangs et al. 2005;Davies et al. 2012), bottom-water warming (Bangs et al. 2005), sedimentation pulse (Zander et al. 2017), and uplift (Foucher et al. 2002). Goto et al. (2016) recently simulated the effect of endothermic cooling from gas hydrate dissociation following depressurisation from tectonic uplift using a temporal thermal sink without modelling gas hydrate dissociation itself. Results predict a significant slowdown of gas hydrate dissociation over thousands of years. Unlike the Goto et al. (2016) study, here we studied the full response of a gas hydrates system using a vertical 1-D TOUGH+HYDRATE v1.5 ("T+H" thereafter) model (Moridis and Pruess 2016;Moridis 2016a) to model gas hydrate dissociation during uplift. This computationally efficient 1-D approach for our conceptual studies allowed for adequate models of fluid pressure changes and fluid flow in the sediment column and testing of a large number of models.

Model Approach
The T+H simulator is an established integral finite difference method simulator developed for the analysis of coupled flow, thermal, and geomechanical processes associated with the formation and dissociation of hydrates in geologic media (Narasimhan and Witherspoon 1976;Moridis and Pruess 2016;Moridis 2016a). It also accounts for all known flow, physical, chemical, and thermodynamic processes associated with the behaviour of hydratebearing systems, where Darcy's law is valid. In this study, our simulated systems are based on a subducting region known to have experienced tectonic uplift, the Nankai Trough offshore Japan. We chose this location because of the availability of relatively well constrained input parameters and for comparison with previous studies (Foucher et al. 2002;Fujii et al. 2016;Goto et al. 2016). Our results, however, are generally applicable to other hydrate settings worldwide locations that are subducting, double or multiple BSR(s) and uplift as been observed, such as the Hikurangi Margin offshore of New Zealand (Pecher et al. 2014(Pecher et al. , 2017.

Model Setup
In our conceptual study, we modelled vertical homogenous 1-D sediment columns and assumed that we initialise for an unconsolidated, unlithified, soft, fully water-saturated marine sediment, and with zero cohesion and tensile strength at a water depth of 600 m, a geothermal gradient of 35 • C/km (Foucher et al. 2002), and at a sea floor temperature of 3.5 • C (Fujii et al. 2016). We also assumed that the initial porosities and intrinsic permeabilities of the sediment columns are constant, like for previous studies (Reagan and Moridis 2008;Liu and Flemings 2009;Thatcher et al. 2013;Marín-Moreno et al. 2015a). In addition, the equilibrium dissociation and formation model was used to account for heat and mass components (water, methane gas, hydrate, and salt) partitioned among four possible phases (gas phase, liquid phase, ice phase and hydrate phase). Using the kinetic or equilibrium hydrate dissociation and formation model should not make much difference because of the time scale of our problem (Moridis and Pruess 2016;Moridis 2016a).
We model the sediment column from the seafloor 0 m (top boundary) to a depth of 1000 m (base boundary) with a 1-D grid of 1,000 elements, using a mesh generator for T+H code called MESHMAKER v1.5 (Moridis 2016b). The elements increase in size with depth, ranging from 0.1 m at the seafloor to 3.711 m at z = 1,000 m. The entire column is equilibrated to initial steady state conditions to ensure stable temperature and hydrostatic pressure distributions and aqueous methane concentrations that correspond to the conditions at the selected depth and temperature (Moridis and Pruess 2016; Moridis 2016a), with a uniform hydrate saturation ( S H ) of 20% within a 20 m-thick hydrate bearing zone above the BHSZ at 160 m, measurement below the seafloor (mbsl), see Fig. 1. Similar input values were used in Goto et al. (2016). 20% is a realistic value for gas-hydrate bearing sands on the Nankai Trough (e.g. Saeki (2008), where saturations of up to 80% were observed (e.g. Matsumoto et al. (2004). We also modelled low-and high-end concentrations of 5% and 50%, respectively.
The top and bottom boundary grid cells of our models were closed, i.e. the boundary conditions were not at constant pressure and temperature. The bottom has a fixed source of heat at 60 mW/m 2 Goto et al. (2016) and fluid flow (saline water) at 1 mm/ yr ( 3.428 ⋅ 10 −7 mol∕m 2 s ) (Pecher et al. 2010), to simulate heat transfer via conduction and advective heat flow typical to subduction zone sediments. To determine the specific enthalpy of the saline water, we setup a simulation with a single grid cell, with initial conditions of prevailing hydrostatic pressure, temperature, and an arbitrary high value for the dissolved methane mass fraction in the saline water. This simulation was performed using a single time step without injection or production, the code printed out the enthalpy for the aqueous phase as 1.715308 × 10 5 J/Kg, mass fraction of dissolved methane and salt in the fluid are at 2.254156 × 10 −3 kg/kg and 0.035 kg/kg. The underlying equations that govern how the code works can be found in Moridis et al. (2019) and Reagan et al. (2019).
The purpose of the dissolved methane mass fraction in the saline water at the base was to overcome the slow diffusive evaporation of hydrate into the overlying seawater (undersaturated with respect to methane) due to the implementation of the diffusion term similar to (Thatcher et al. 2013;Marín-Moreno et al. 2015. The intrinsic permeabilities used for this study are 10 −13 , 10 −15 , and 10 −17 m 2 simulating sands, silts, and clays, respectively (Spinelli et al. 2004;Matsumoto et al. 2004). Gas saturation (S g ) below BHSZ is fixed at 2 %, consistent with previous studies (Marín-Moreno et al. 2015a;Thatcher et al. 2013;Reagan and Moridis 2008). We used the modified version of Stone's first three-phase relative permeability model to calculate the relative permeability for methane and water (Stone et al. 1970). Capillary pressure was computed using the van Genuchten function (Van Genuchten 1980) and for changes in capillary pressure and intrinsic permeability due to changes in porosity, pressure and hydrate saturation in pore spaces, we used the evolving porous media method (EPM #2) (Moridis and Pruess 2016; Moridis 2016a).

Top Boundary Condition
In our T+H model, the water column is represented as the pressure at a point on the first grid cell and as a boundary and is treated as an "infinite volume". Boundary conditions are unaffected by what happens within the domain: they affect it but are not affected by it (this is the concept of a boundary). We used a time-variable boundary, i.e. the time-dependence is read from a table. There is no attempt to re-calibrate the pressure because that would be nonphysical: therefore, if the boundary pressure changes, then the pressure in the entire domain changes (Moridis and Pruess 2016; Moridis 2016a). Braga et al. (2022) employed this approach to evaluate the response of shallow methane hydrates to simultaneous sea level increases during the last glacial maximum on the Amazon Deep-Sea Fan in Brazil.
For the actual simulation, the seafloor boundary changes when the pressure at the seafloor changes to reflect changing sea level, the distance between the water surface and the seafloor. Hydrostatic pressures are then calculated for the remaining grids representing the marine sediment. We modelled the pressure change at the seafloor boundary by allowing it to propagate very fast and dynamically depending on the simulation time scale, since change in pressure (sea level) does not occur instantaneously in nature (Hermanrud et al. 2013). However, there is a small but finite time for the pressure change to propagate into the sediments based on Darcy flow, which is similar to other numerical studies, e.g. (Bredehoeft and Papaopulos 1965;Bredehoeft and Hanshaw 1968;Hermanrud et al. 2013). Other physical properties factored into how the pressure change propagates into the sediment column include: porosity, permeability, fluid gravity, pressure, temperature, geomechanical effects, and thermal properties of the sediments in relation to any other multiphase behaviour (Moridis and Pruess 2016; Moridis 2016a).

Modelling Uplift
In our study, we simulated the uplift by changing the fluid pressure of the whole sediment column by decreasing the pressure at the seafloor (the top of the simulated sediment column), which is a similar approach used in the Goto et al. (2016) tectonic uplift study. We decreased pressure at the seafloor at a rate of 0.009 myr −1 over a duration of 8kyrs (8,000 years), at a constant seafloor temperature based on data from Foucher et al. (2002). We consider 8 kyrs as the key time frame for investigating uplift because 7-8 kyrs ago, the postglacial eustatic sea-level rise decelerated significantly (Lambeck et al. 2014). Before then, the increase of pressure from the sea-level rise at rates of >14 mm/yr (Fleming et al. 1998) more than counteracted any depressurisation from uplift. This approach has also been taken by (Pecher et al. 2017) to deduct uplift rates from d-BSRs on the Hikurangi Margin.

Poro-Elastic Effect
We also took note of the poro-elastic effect of gas-bearing sediments on pore pressure using the loading efficiency ( LE ), which controls how much change in fluid pressure is caused by an increase/decrease in the applied load under undrained conditions assuming an incompressible sediment solid (Wang 2000;Sawyer et al. 2008;Tanikawa et al. 2008;Liu and Flemings 2009;Hermanrud et al. 2013). The magnitude of the pore pressure response depends on the compressibilities of the fluid and porous matrix and is applicable to a sudden or rapid change in the mean stress (Wang et al. 1998;Wang 2000;Katahara and Corrigan 2001;Sawyer et al. 2008;Liu and Flemings 2009).
According to Sawyer et al. (2008), a sediment containing no free gas as a loading efficiency between 0.95 and 0.99, implying an increase/decrease in pore pressure equal to an increase/ decrease in the applied load with a slight change in the effective stress. This effect varies over time and space depending in particular on gas saturation, since gas is highly compressible. The T+H code cannot account for this effect during its model run, which is a limitation of the code. Our modelling stretches over longer time spans which should allow the models to reach equilibrium and thus, we do not expect the poro-elastic effect to affect our results significantly. To test this assumption, we ran a model simulating the poro-elastic effect of constant S g at 2%. We also assumed a loading efficiency of 0.50, this is because of the presence of free gas in the sediment, which increases the pore fluid compressibility and lowers loading efficiency between 0.20 -0.80 (Wang et al. 1998;Wang 2000;Liu and Flemings 2009).
According to Liu and Flemings (2009), LE can be described using the equation below: where m v is the formation compressibility (where m v = −dV∕Vd , V is the sediment volume and is the vertical effective stress), is the porosity, S w is water saturation, w is the water compressibility, g is the gas compressibility and is described by: In Eq. 2, P is the pressure and Z is the gas compressibility factor. For an ideal gas case, Z = 1 and g is the reciprocal of the pressure. In addition, to account for the change in pore pressure due to poro-elastic effect in our study. We determined the formation compressibility of the sediment using 4.5 × 10 −10 Pa −1 as w (Yamamoto et al. 2017), S w is 98%, and g at 600 water depth.

Pore Pressure Parameter
The dissociation of gas hydrate into free gas and water in pore spaces of a rock matrix may produce elevated pore fluid pressure (Flemings et al. 2003). However, the rate of elevated pore fluid pressure depends on the rate of gas hydrate dissociation and the permeability of the sediment (Xu and Germanovich 2006). Consequently, an increase in pore fluid pressure could either increase the fluid mobility or fracture the rock matrix (Daigle and Dugan 2010). To analyse the elevated pore-pressure, we computed the pore pressure parameter ( * ) using the equation by Screaton et al. (1990): where P is the computed pore fluid pressure, P hydro is hydrostatic fluid pressure and P lith is the lithostatic pressure. If * is equal to zero, pore pressure is hydrostatic, if it is one, it is lithostatic. * ≥ 0.6 is often used for the occurrence of vertical hydrofracturing (e.g. (Daigle and Dugan 2010)), i.e. when the pore pressure in the rock matrix is greater than the horizontal confining pressure, see Table 1 for more parameters and constrained used in this study.

Results
We simulated the response of gas hydrates in sediments to 0.009 myr −1 rate of tectonic uplift (Foucher et al. 2002). Figure 2 shows the hydrate, gas, and water saturation profiles for an initial S H = 20%, while Fig. 3 shows the hydrate profiles for an initial S H = 5% and 50% during uplift. Figure 4 shows the pore pressure distribution profiles for an initial S H = 5%, 20 %, and 50 % during uplift, while Fig. 5 shows the salt mass fraction, pore pressure parameter, and temperature deviation during uplift for an initial S H = 20%. We first evaluate the sands and silts scenarios (see Fig. 2) since the general pattern of hydrate, gas, and water distribution is similar (see Fig. 2a and b, d and e, and g and h). Two Table 1 Physical rock properties and simulation parameters a S A , S G , and S H are saturations for aqueous, gas, and hydrate phases. S irG and S irA are irreducible gas and aqueous saturations and S mxA is the maximum water saturation. K rG and K rA are relative permeabilities for gas and aqueous phases. P cap is the capillary pressure, P 0 is the maximum value of capillary pressure, and P max is the capillary entry pressure; and n, and are fitting parameters 1 3 separate hydrate layers are developing, a layer of dissociating gas hydrate within the original depth interval of hydrate occurrence as well as a newly forming hydrate layers above that, see Fig. 2a and b. Figure 3 shows the profiles of gas hydrate saturation for an initial S H of 5% and 50%. The 50% hydrate saturation model's general pattern for hydrate distribution is similar to that of 20% initial S H . For the 5% case, the lower hydrate layer for the sands and silts scenarios disappeared after 6kyrs and only one layer is observed after 8kyrs.
The gas saturation S G is zero within the zone of gas hydrate occurrence for the initial model before uplift. During uplift, gas forms within the original layer of hydrate occurrence and moves gradually upwards (Fig. 2d, e). The position of the gas front coincides with the top of hydrate occurrence and thus, gas is present in the entire depth interval across the base of the original gas hydrate layer to the top of the newly forming hydrates. Water saturation (Fig. 2g, h) mimics the evolution of gas hydrate dissociation and formation. Figure 4 shows the pore pressure distribution profiles for sands, silts, and clays scenarios during uplift. We observe that for the duration of our sands (see Fig. 4a, d, and g), silts (see Fig. 4b, e, and h) simulations, the pore pressure remained hydrostatic for all hydrate saturations (5%, 20%, and 50%). The pore pressure parameter for sands and silts is not elevated significantly (see Fig. 5d, e).
The salinity profiles are used as a proxy for gas hydrate formation or dissociation: ion exclusion during hydrate formation increases salinity in the remaining pore water, while the release of fresh water from dissociating hydrate decreases salinity. The sands a b c d e f g h i Fig. 2 Evolving hydrate a, b and c, gas d, e, f and water g, h and i saturation profiles for sands, silts and clays for an initial hydrate saturation of 20% with basal fluid flow. The z-axis for each graph represents depth below the seafloor. The black lines represent the initial state of each observed saturation and silts scenarios (see Fig. (5a, b) show two types of excursions from across the black curve which represent the initial molecular weight of the dissolved salt mass fraction (0.035 wt %) (Reagan and Moridis 2008). A decrease in salinity is evident in the original zone of gas hydrate occurrence whereas salinity increases in the zone of newly formed hydrate. Temperature is displayed as deviation from the initial temperature gradient of 35 • C/km, see Fig. (5g, h, and i). For the low-permeability clays scenario, an entirely different picture emerges (see Fig. 2c, f, and i). Gas hydrate remains in place at its original position after 8 kyrs, dissociating to gas and water. The gas does not appear to migrate upwards. An elevated pore pressure was observed below 140 m (i.e. from the top of the hydrate layer) for all hydrate saturations (5%, 20%, and 50%) (see Fig. 4c, f, and i). For the simulation containing hydrate saturation of 50%, the dissociation of gas hydrates leads to an increase of pore pressure above the original hydrostatic pressure below 160 m (see Fig. 4i). This is caused by the presence of gas released from hydrate dissociation during uplift. The pore pressure parameter increases from 0 -∼0.40 with uplift indicating pore pressure buildup (Fig. 5f). The salinity profile shows two major layers of slight decrease in salinity (see Fig. 5c), indicating very slow dissociation of hydrates with most hydrate remains stable (see Fig. 2c) through the duration of uplift. The zones of slightly elevated salinity beneath the BHSZ and top of hydrates (see Fig. 5c), confirms a formation of hydrate beneath the BHSZ (see Fig. 2c). This is due to the pore pressure increase caused by the advective fluid flow at the base of our model (see Fig. 4f). The temperature drop around the location of hydrate for both the sands and silts scenarios, see

Sensitivity Tests for Poro-Elasticity and Basal Fluid Flow
The models above did not account for any poro-elastic effects. While T+H cannot account for poro-elastic during the model run, we investigated its potential impact by assuming a loading efficiency of 0.50. Figure 6 shows the hydrate and gas saturations, and pore pressure distribution profiles for sands, silts, and clays scenarios for the response of poro-elastic effect on our model (Liu and Flemings 2009). The hydrate saturation profiles for sands, silts, and clays in Fig. 6a, b and c show similar patterns to Fig. 2a, b and c. Similar partners are observed for the gas saturation profiles, see Fig. 6d, e and f show similar patterns to Fig. 2d, e and f. For the clays scenario, the elevated pore pressure seems to be lower than for the case where poro-elasticity has been not been considered, see Fig. 6i and Fig. 4f, while the sands and silts cases show similar patterns where poro-elasticity has been not been considered, see Fig. 6g and h and Fig. 4d and f. We further injected the dissolved methane mass fraction only, without additional fluid flow, to investigate the role of basal advective fluid flow in our model like in previous studies (Liu and Flemings 2009). Therefore, heat transfer is via conduction only. For the sands and silts scenarios (see Fig. 7), the distribution of hydrate, in Fig. 7a and b, gas in Fig. 7d and e, and water in Fig. 7g Fig. 2a and b, gas in Fig. 2d and e, and water in Fig. 2g and h saturations. Similar behaviours are observed for the salinity in Fig. 8a and b, pore pressure parameter in Fig. 8d and e, temperature deviation in Fig. 8g and h, and pore pressure profiles in Fig. 9a, b, d, e, g and h and for models with basal advective fluid flow for salinity in Fig. 5a and b, pore pressure parameter in Fig. 5d and e, temperature deviation in Fig. 5g and h, and pore pressure profiles in Fig. 4a, b, d, e, g and h profiles.
For the lower-permeability clays scenario, a different behaviour was observed compared to the sands and silts scenarios. Gas hydrate continued to dissociate from the base of hydrate stability for the duration of 8 kyrs, see Fig. 7c, leading to the release of gas, see Fig. 7f and water, see Fig. 7i. On the other hand, hydrates were observed to be in place for 8 kyrs for the model with basal advective fluid flow, see Fig. 2c. The gas released behaved in a similar manner as the model with basal advective fluid flow, see Fig. 2f, i.e. the gas does not migrate upwards. However, the gas released from hydrate is ∼12% more than the initial saturation, compared between Fig. 7f and Fig. 2f. The fresh water released, see Fig. 7i, from gas hydrate dissociation causes salinity to decrease, see Fig. 8(c), further than what is observed in the model with basal advective fluid flow, see Fig. 5c. Furthermore, the pore pressure is not as elevated as that of the model with advective fluid flow, see Fig. 9f and see Fig. 5f for comparison.  salinity (a, b and c), pore pressure parameter (d, e and f), and temperature deviation (g, h and i) profiles for sands, silts and clays for an initial hydrate saturation of 20% with basal fluid flow. The z-axis for each graph represents depth below the seafloor. The black curve represents the initial state of each observed quantity. Please note the different scales for the salt mass fraction for the clays scenario

Discussion
Our key finding is that gas hydrate may remain stable close to its pre-uplift BHSZ for thousands of years. The dissociation of methane hydrate is an endothermic reaction, that is, it absorbs heat from nearby sediment to dissociate gas hydrate into water and gas. This process causes the local temperature to decrease at the BHSZ and thus, stabilise the gas hydrate. Conductive and advective heat transport is too slow to dissipate the resulting thermal anomalies entirely over 8 kyrs. This result is similar to the findings reported in thermal models that treated gas hydrate dissociation as a heat sink and that the latent heat of gas hydrate dissociation plays a significant role in the process of dissociation of gas hydrate (Goto et al. 2016). A decrease in salinity and increase in pore pressure from gas hydrate dissociation further stabilise gas hydrates resulting in slowing down gas hydrate dissociation.
For higher permeabilities of the sands and silts scenarios, we assumed that pore spaces are well connected within the sediment. Therefore, the elevated pore pressure caused by gas hydrate dissociation and basal advective fluid flow is every small (Xu and Germanovich 2006;Goto et al. 2016) and dissipates very quickly, see Figs. (5d and e and 8d and e). Gas and water generated from gas hydrate dissociation migrate upwards, see Figs. (2d, e, g and h and 7d, e, g and h). We also observed that the percentage of initial hydrate dissociation decreases with increasing hydrate saturation, see Evolving hydrate a, b and c, gas d, e and f, and pore pressure distribution g, h and i profiles for sands, silts and clays for an initial hydrate saturation of 20% with basal fluid flow and the poro-elastic effect for a constant S g at 2%. The z-axis for each graph represents depth below the seafloor. The black lines represent the initial state of each observed saturation and b and 3a, b, d and e). However, we observed that the combination of advective and conductive heat flow causes pre-uplift hydrates to dissociate faster, see Fig. 2a and b than in models with only conductive heat flow, see Fig. 7a and b. A second shallower layer of hydrate forms in the sediment column. This second hydrate layer is thicker in models with advective and conductive heat flow, see Fig. 2a and b, while models with only conductive heat flow have higher hydrate saturation , see Fig. 7a, b.
Ongoing hydrate dissociation and reformation are reflected by a decrease and increase of salinity and temperature deviation, respectively, see (Figs. 5a, b, g and h and 8a, b, g and h). The initial temperature gradient at the seafloor is slightly elevated, which we attribute to the contribution from the heat source below and possibly advective and conductive heat flow from that the heat source leading to an increase of advective heat flux due to the additional fluid source from hydrate dissociation, see Fig. 5g and h. Alternatively, the conductive heat flux may increase temporarily because exothermic hydrate formation as a heat source occurs closer to the seafloor than endothermic hydrate dissociation. In comparison, the model without advective heat flow, we consistently see a decrease of the geothermal gradient at the seafloor, see Fig. 8g  The modelled spikes in hydrate saturation in the sands scenario may not be an artefact but could mark locations where local pressure-temperature conditions cross the phase boundary, see . (2a and 7a). This could be caused by advective fluid flow and quick formation of hydrates, which leads to an increase in salinity (see (Fig. 5a) and temperature within the hydrate stability zone (see (Fig. 5g). However, the original in situ pore pressure is less than the pore pressure conditions needed for hydrate stability, making the hydrate less stable and it forms elsewhere (see (Fig. 4d). Hydrate formation is sufficiently slow that free gas, water, and hydrate co-exist. This could be due to the modelled increase in salinity as suggested by other models Flemings 2006, 2007) and field observations Torres et al. 2004). The upper level of hydrate occurrence is determined by the top of the gas front (Fig. 3d, e).
For the low permeability of the clays scenario, we assumed that our sediment's pore spaces are confined. Therefore, the elevated pore pressure caused by basal advective fluid flow and gas hydrate dissociation causes an increase in volume expansion and does not dissipate quickly, leading to pore pressure build-up (Xu and Germanovich 2006;Goto et al. 2016), see Fig. 5f. During repeated uplifts, the gas hydrate barely dissociated. Any gas released from hydrate due to uplift does not migrate upwards, see Fig. 2f. Instead, the pore pressure around the gas hydrate bearing region is elevated to the initial hydrostatic pore pressure, see Fig. 4f, and for higher S H , e.g. 50%, the pore pressure was elevated above the initial pore pressure, see Fig. 4i. On the other hand, the models without basal advective fluid flow show a slight pore pressure increase (see Fig. 9f i profiles for sands, silts and clays for an initial hydrate saturation of 20% without basal fluid flow. The z-axis for each graph represents depth below the seafloor. The black curve represents the initial state of each observed quantity. Please note the different scales for the salt mass fraction for the clays scenario to gas hydrate dissociation (see Fig. 7c) leading to pore pressure build-up (see Fig. 8f).
In addition to the effect of endothermic cooling caused by uplift (Goto et al. 2016), hydrate dissociation is thereby slowed down, see Fig. 5i. The low permeability furthermore slows down the diffusion or advection of locally decreased salinity from hydrate dissociation further stabilising gas hydrates, see Fig. 5c. The pore pressure parameter is noticeably elevated, although still below the "rule-of-thumb" value of 0.6 for hydrofracturing (Daigle and Dugan 2010) (see Fig. 5f). This study has significant implications for gas hydrate bearing systems in tectonically active regions. Most importantly, we predict a layer of dissociating gas hydrate to be present for thousands of years at the level of the predicted pre-uplift BHSZ at the end of the postglacial sea-level rise (see Fig. 10). Endothermic cooling from gas hydrate dissociation is predicted to lead to a depressed geothermal gradient even for seafloor thermal measurements. For the higher permeabilities in silts and sands, gas and hydrate coexist over much of the depth range from the original base of hydrate occurrence to the top of newly forming gas hydrates. Elevated pore pressure from hydrate dissociation dissipates, for the low-permeability clays scenario, the gas remains in place after dissociation. Our models show an increase in pore pressure for models with basal advective heat flow, (see Fig. 4c, f, and i), while Fig. 9c, f, and i shows similar models without basal fluid flow. However, if the elevated pore pressure already existed before uplift, it is possible that the additional pore pressure from hydrate dissociation and advective heat flow leads to a decrease of sediment strength. Evolving pore pressure distribution a, b and c, d, e and f, and g, h and i profiles for sands, silts and clays for an initial hydrate saturation of 5 %, 20%, and 50% without basal fluid flow. The z-axis for each graph represents depth below the seafloor. The black lines represent the initial state of each observed saturation The introduction of the poro-elasticity effect to our initial S H = 20 % model (see Fig. 2) is unlikely to change anything significantly. Similarly, advective flow does not have any major effect on our models, see Fig. 6. However, the pre-uplift hydrates dissociate faster in the sands and silts scenarios compared to our original model, see Fig. 6a and b and 2a and b. The clays scenario shows hydrate dissociating from the base of hydrates compared to our original model, see Fig. 6c and Fig. 2c for comparison. The thickness of methane gas released from hydrate dissociation is slightly higher in the sands and silts scenarios compared to our original model, see Fig. 6d and e and Fig. 2d and e, while the clays scenario shows more gas released than in the original model, see Fig. 6f and Fig. 2f. Lastly, the pore pressure profiles show the same pattern, implying we would be over-estimating or underestimating the pore pressure but will arrive at similar results, indicating this limitation to have a negligible effect on our model, see Fig. 6g, h and i, Fig. 4d, e, and f, and Fig. 9d, e, and f.
Our results have implications for the interpretation and nature of BSRs and double-BSRs in regions of tectonic uplift. For the sands and silts scenarios, the two hydrate layers would intuitively be compatible with the presence of d-BSRs. However, the presence of co-existing gas and hydrate makes it difficult to predict the seismic response. Fig. 10 Schematic sketch not to scale showing the effect of uplift on high-and low-permeable marine hydrate bearing sediment for a duration of 8 kyrs. In the high-permeable marine hydrate bearing sediment two separate hydrate layers were observed in a gas saturated sediment. While in the low-permeable marine hydrate bearing sediment only one hydrate layer is observed Most rock physics models for gas hydrates in nature predict the drop of P-wave velocity ( V p ) from gas even at a saturation of only a few percent to be much more significant than the V p increase from gas hydrate at saturations below ∼ 40% (e.g. (Helgerud et al. 1999;Bai et al. 2016))-the only exception would be hydrate cementation of grain or between grain contacts, which has to our knowledge only been documented at low saturations in laboratory experiments (e.g. Priest et al. (2005)). For the typical values of gas and hydrate saturations in our models, we would therefore expect a V p decrease and thus, a negative polarity in seismic reflections similar to that from BSRs, at the top of co-existing hydrate and gas occurrence. We speculate that in the seismic data, the interval of coexisting gas and hydrate could be interpreted as a free gas zone beneath the BHSZ. Internal reflectivity within this zone of coexisting gas and hydrate may be linked to changes in hydrate saturation.
If the permeability is too low such that gas cannot migrate (our clays scenario), a negative polarity reflection will be expected to remain in place at the paleo-BHSZ. Our models predict that no new gas hydrate will form above the original zone of hydrate occurrence. In reality, we speculate that a new BSR may form at the post-uplift level of the BHSZ from already existing gas hydrates, perhaps in different lithologies that allow gas to migrate or perhaps by gas supplied along the fractures rather than through pore-scale permeability. In the low-permeability case, our models largely confirm the presence of a depressed level of one BSR. However, previous concepts of uplift-related paleo-BSRs invoke the gas itself to be trapped by low relative permeability at the level of the paleo-BHSZ while the hydrate dissociates (Bangs et al. 2005). We now predict gas hydrate to remain in place at its pre-uplift level, not just gas.

Conclusions
A drop in pressure caused by tectonic uplift leads to the dissociation of methane hydrate into water and gas. For our more permeable scenarios mimicking sands and silts, elevated pore pressure dissipates, fluids migrate towards the seafloor and a shallower hydrate layer forms. For the low-permeability clays scenario, the gas released from hydrate remains within the original zone of hydrate occurrence. Endothermic cooling, along with salinity decrease and elevated pore pressure, slows down hydrate dissociation to a degree that dissociating gas hydrate exists for thousands of years beneath the post-uplift regional level of the BHSZ. While our 1-D models are not suitable for modelling the formation of the shallower BSR, we predict the deeper BSR to mark gas that remains in place beneath slowly dissociating hydrate at the pre-uplift level of the BHSZ.