Microbial Methane Generation and Implications for Stability of Shallow Sediments on the Upper Slope, U.S. Atlantic Margin

Dissociation of methane hydrates in shallow marine sediments due to increasing global temperatures can lead to the venting of methane gas or seafloor destabilization. Along the U.S. Atlantic margin there is a well-documented history of slope failure and numerous gas seeps have been recorded. However, it is not fully understood whether the observed gas seepages can lead to slope failure as seafloor data is often sparse. We used machine learning algorithms to predict total organic carbon (TOC) and porosity at the seafloor on the U.S. Atlantic margin. Within this region, an area of high TOC predictions (1.5—2.2% dry weight) occurred along the continental slope from (35.4°N, 75.0°W) to (39.0°N, 72.0°W), aligning with documented gas seeps in the region. Elsewhere, predicted values of TOC were near or below 1% dry weight. In the area of high TOC, we modeled hydrate and gas formation over a 120,000 years glacial cycle. Along the feather edge, average hydrate saturations at the base of the hydrate stability zone (BHZ) were between 0.2% and 0.7% with some models predicting hydrate saturation above 3% and average peak gas saturations ranged from 4% to 6.5%. At these locations we modeled the pore pressure response of sediments at the BHZ to hydrate dissociation due to an increase in temperature. We focused on purely drained and undrained loading environments and used a non-linear Hoek-Brown failure envelope to assess whether failure criteria were met. In a drained loading environment, where excess pore pressure is instantly dissipated, we found that the change in effective stress due to hydrate dissociation is small and no failure is expected to occur. In an undrained loading environment, where excess pore pressure does not dissipate, the change in effective stress due to hydrate dissociation is larger and shear failure is expected to occur even at low hydrate saturations (0.2%—1%) forming final gas saturations below 0.1%. Therefore, we conclude that the dissociation of hydrates along the feather edge can lead to the conditions necessary for sediment failure.


INTRODUCTION
Methane hydrates are solid, crystalline structures composed of methane and water. Hydrates are stable in an environment with low temperature, high pressure, and where sufficient methane is present (Kvenvolden and Claypool, 1988). Hydrates are stable beneath the seafloor when the bottom water temperature approaches 0°C and the water depth is greater than 300 m (Kvenvolden and Lorenson, 2001). The depth range at which methane hydrates are stable is known as the gas hydrate stability zone (GHSZ). The GHSZ continues below the seafloor to the base of the hydrate stability zone (BHZ) where the hydrate phase envelope intersects the geotherm. The depth of the BHZ can be found up to 2000 m below the seafloor but is typically found at much shallower depths (Kvenvolden, 1993).
An increase in local temperature or decrease in local pressure can cause destabilization and the subsequent dissociation of hydrates (Phrampus and Hornbach, 2012;Ruppel and Kessler, 2017). Methane gas released due to hydrate dissociation can migrate up to the seafloor and into the water column. In the water column, it can oxidize into CO 2 causing ocean acidification, or it may continue through the water column and into the atmosphere as a greenhouse gas (Biastoch et al., 2011). Increasing temperatures and the subsequent dissociation of hydrate is especially a concern at the BHZ where the in-situ pressure and temperature conditions intersect the hydrate phase envelope and even minor changes in temperature can cause hydrate dissociation. There is evidence that within the last 100 years, increases in global ocean temperatures have led to hydrate instability and the shifting of the hydrate stability zone along the Beaufort shelf and off the Svalbard coast (Ferré et al., 2012;Sarkar et al., 2012;Phrampus et al., 2014). In shallower water (around 500 m) where the BHZ is at or near the seafloor, also known as the feather edge, the increase in seafloor temperature is problematic (Ruppel, 2011). As water depth increases, this becomes less of an issue as the GHSZ thickens. In shallow marine sediments within 132 m of the seafloor, the presence of gas can lead to tensile fracturing in the sediments, and gas in these shallower sediments can occupy a much larger volume than hydrate under the same conditions (Daigle et al., 2020).
We are interested in the geomechanical properties of the seafloor sediment when methane gas is formed due to hydrate dissociation from increasing temperatures. We focus specifically on the feather edge where the hydrate stability zone is the thinnest. To begin, we model TOC along the U.S. Atlantic margin. At a few locations where estimated seafloor TOC is high, we use a 1D sediment burial and methanogenesis model to simulate hydrate and gas formation over a 120,000 years glacial cycle. At locations along the feather edge, we assume a 1°C increase in temperature at the BHZ and model the geomechanical response of the system as hydrate dissociates to methane gas and water. We are interested in determining if the conditions for sediment failure due to an increase in temperature can exist under the right circumstances. Thus, at each location we calculate a Hoek-Brown failure envelope to determine if failure occurs due to dissociation.

BACKGROUND
Ocean Drilling Program (ODP) boreholes give an insight into locations where methane gas and methane hydrate occur along the U.S. Atlantic margin. Along the upper continental slope, offshore of New Jersey, ODP boreholes report subsurface microbial methane, indicating methanogenesis in the area (Shipboard Scientific Party, 1994a;Shipboard Scientific Party, 1994b;Shipboard Scientific Party, 1994c;Shipboard Scientific Party, 1994d;Shipboard Scientific Party, 1994e;Shipboard Scientific Party, 1998a). South of this region, offshore of North Carolina and South Carolina, seismic profiles from the ODP show a strong bottom-simulating reflector (BSR) at some locations, inferring the presence of hydrate in the area (Shipboard Scientific Party, 1996b;Shipboard Scientific Party, 1996c;Holbrook et al., 1996;Dickens et al., 1997).
In sediment near the seafloor, gas-driven tensile failure can occur due to a combination of low stresses and weak sediments. Fine grained sediments with a higher clay fraction are more prone to failure as gas within these sediments has a higher capillary pressure. Therefore, less gas is required to be present to cause failure (Daigle et al., 2020). One mechanism for gas generation is the dissociation of methane hydrates although gas generated directly from microbial methanogenesis can also lead to tensile failure.
The release of methane gas from the seafloor can be seen in water column acoustic anomalies which correspond to gas plumes (Judd, 2003;Skarke et al., 2014). In global surveys summarizing free gas distribution in marine sediments, Fleischer et al. (2001) and Judd (2003) both note occurrences of free gas along the U.S. Atlantic margin. The release of gas is evident in the development of pockmarks on the seafloor which form due to the accumulation of gas beneath a seal, and the subsequent release of gas when the seal fails (Cathles et al., 2010;Sultan et al., 2014). Along the U.S. Atlantic margin, pockmarks just shallow of the feather edge were reported by . In addition, asymmetric depressions corresponding to elongated gas blowouts were reported by Hill et al. (2004) along the U.S. Atlantic margin.
Investigations by Skarke et al. (2014) found instances of methane gas leakage from the seafloor along the U.S. Atlantic margin at a higher concentration than previously thought. Skarke et al. (2014) identified 570 gas plumes on the northern U.S. Atlantic margin, with 440 of these seeps (77%) lying between the shelf break (~180 m below sea level) and 600 m below sea level. The location of these plumes would lie just updip of the feather edge of the GHSZ. The seeps in this area were further explored by Prouty et al. (2016) who suggest that seepage may have begun as early as 15 kya in the Baltimore Canyon seep field, and between one to three kya at the deeper Norfolk seep field ( Figure 1). Gas seeps in this region have δ 13 C values between -73.5‰ and -109‰ and methane:ethane ratios between 385 and 926 indicating that the gas is microbial in origin (Pohlman et al., 2017).
Dakota and PFLOTRAN to model the generation of methane gas and hydrate due to methanogenesis and the burial of the seafloor sediment. In locations where the BHZ was close to the seafloor, we modeled the geomechanical response of sediment to gas generated by hydrate dissociation to better understand whether shear failure or tensile failure would occur. We specifically focused on hydrate dissociation caused by an increase in temperature at the BHZ. The hydrate formed through the PFLOTRAN model and subsequently dissociated in the geomechanical model was assumed to be structure I methane hydrate.

The Global Predictive Seabed Model
GPSM is a geospatial machine learning model developed by the Naval Research Laboratory and can be used to predict TOC, porosity, sediment composition, and other properties of seafloor sediments in areas where no measurements have been sampled or recorded (Martin et al., 2015). GPSM offers a variety of machine learning methods to create geospatial models including k-nearest neighbor regression, support vector regression, and random forest regression. In these geospatial machine learning models, predictions do not assume spatial autocorrelation. Instead, the characteristics or parameters of locations are compared rather than the geographic proximity of two points.

Modeling Seafloor Total Organic Carbon
To model seafloor TOC, we followed the methodology of Lee et al. (2019) and used a k-nearest neighbor (kNN) approach. For our model, a value of k = 5 was chosen. When predicting TOC values on the seafloor at unknown locations, GPSM compares various observed attributes of the unknown location with the same observed attributes at locations with known TOC values. These observed attributes are known as predictors and are found in global grids known as predictor grids. In our model 1749 predictor grids were used including data from multiple global surveys such as seamount censuses, global river fluxes of carbon and sediments to the ocean, and decadal trends in oxygen concentration in subsurface waters, among others . Therefore, our model differs slightly from that of Lee et al. (2019) as we are using both new and additional predictor grids.
In our model, 5,595 individual locations with known TOC values were sampled (Seiter et al., 2004). These are from the upper 5 cm of sediment and so are roughly Pleistocene in age. Of these, we excluded values of TOC over 5% to mitigate outliers in the predictions. We then averaged values for each 5 × 5 arc-minute grid cell to obtain a uniformly spaced grid. Of the 5,595 locations with known TOC, 126 points had measurements greater than 5% TOC (2.2% of the data) and were excluded. After creating a uniformly spaced grid, 4,879 useable observations remained ( Figure 1).
We used tenfold validation to validate our model. In tenfold validation, the data is first randomly split into 10 equal sized groups. One group is excluded as a test set and the remaining 90% of the data is used to form a training set. A model (in this case kNN) is created using the training set of data and then tested on the remaining 10% of the data where the predicted value from the model can be compared to the actual value of the test set. This process is done 10 times. Plotting the observed versus predicted data values gives us the validation graph in Figure 2. The validation of this data shows an R 2 value of 0.6151, lying below the ideal 1:1 fit at high observed TOC values. Other k values between 3 and 20 were tested to optimize for the highest R 2 value. We chose k = 5 for our model as it was among the top performing models and is consistent the work done by Lee et al. (2019). After forming a model using tenfold cross validation, we checked the predictor grids for collinearity. Grids with a correlation coefficient over 0.99 were removed, leaving 1,622 of the 1,749 predictor grids. In addition, an error grid was used to remove predictor grids with high error. Since this grid is made of random, uniform noise, it should have no influence on the data being modelled. The predictors grids can be ordered by individual error and any predictor grids with a higher error than the noise grid is ignored and not used in the model. This left us with 1,618 predictor grids which were used for the prediction. TOC predictions were made on the 5 × 5 arc-minute grid within the area 29°N-45°N and 82°W-66°W ( Figure 1).

Modeling Seafloor Porosity
We also used GPSM to model seafloor porosity in our region of interest. The same 1,749 predictor grids used to predict seafloor TOC were used in our seafloor porosity model. In this model, 2027 individual locations with known seafloor porosity values were sampled (Martin et al., 2015). Of these, values of porosity over 90% and below 20% were excluded and the remaining values were averaged for each 5 × 5 arc-minute grid leaving us with 1,440 useable observations. Both a kNN approach and a random forest regression (RFR) approach were modeled. The RFR methodology followed a similar workflow to the work done by Graw et al. (2021) and is similar to the work done to predict seafloor TOC. As we previously did, we used a tenfold cross validation method to determine model variation. A similar validation plot to TOC was calculated for both the kNN and RFR methods using tenfold validation for seafloor porosity. A maximum R 2 value for the RFR method was calculate using 100 trees and a minimum samples per leaf of 3 ( Figure 2). Although this R 2 value (0.4185) was less than the R 2 value for the kNN test (0.4562), we chose to use the RFR model due to lower standard deviations of the porosity predictions at our points of interest.

Dakota and PFLOTRAN
Dakota and PFLOTRAN were used to model the burial of seafloor sediments and the resulting generation of methane hydrate and methane gas in an approach similar to Eymold et al. (2021). Both Dakota and PFLOTRAN are open-source software developed by Sandia National Laboratories. Dakota was developed to provide optimization tools for simulations and we used it to sample PFLOTRAN input parameters, creating a distribution of model results (Adams et al., 2021).
PFLOTRAN is a parallel subsurface flow code which solves a system of nonlinear partial differential equations describing flow and transport in porous medium (Hammond et al., 2014). Following the work of Eymold et al. (2021) and Nole et al. (2017), we used PFLOTRAN to simulate hydrate and gas formation from microbial gas during sediment burial. Eymold et al. (2021) and Nole et al. (2017) use a primary variable switch method to solve for mass and energy conservation with phase change (hydrate versus gas phase). The governing equations are set up as a finite volume difference discretization, and variables are then solved for in a fully implicit manner using a nonlinear Newton-Raphson iteration. The methane generation rate was calculated following Malinverno (2010): where q(z) is methane generation rate with depth in kg/m 3 /s, k α = 2,241 kg/m 3 is the conversion factor of TOC to methane, λ is the methanogenesis rate (s −1 ), ω is the sedimentation rate (m/s), and z is the sediment depth (m). The depth of the sulfate-methane transition (z SMT ) was set at 15 m below seafloor (mbsf) based on work done by Malinverno (2010) and Egger et al. (2018). The metabolizable organic carbon remaining at the SMT (α) was set at 75% (Bhatnagar et al., 2007). Initial conditions and boundary conditions for the simulation were also set up following the methodology of Eymold et al. (2021) and Nole et al. (2017). The initial conditions were set using Dirichlet temperature values, Dirichlet mole fraction values, and hydrostatic pressure at the seafloor. Temperature was set throughout the sediment column at a fixed temperature gradient from the seafloor. In addition, the mole fraction for methane in the aqueous phase was set to 0.001 along the entire profile. At the U.S Atlantic margin, there is documentation of gas in sediments as old as the Oligocene as well as deeper Mesozoic accumulations of thermogenic gas (Shipboard Scientific Party, 1994a; Shipboard Frontiers in Earth Science | www.frontiersin.org May 2022 | Volume 10 | Article 835685 Scientific Party, 1994b;Shipboard Scientific Party, 1994c;Shipboard Scientific Party, 1994d;Shipboard Scientific Party, 1994e;Party, 1998a;Shipboard Scientific Party, 1998a;Party, 1998b). However, since we were most interested in gas and hydrate close to the seafloor and the role of seafloor TOC in generating it, we chose to ignore the hydrate and gas that may already exist in the areas we model. At the bottom of the depth profile, Neumann boundary conditions were set such that the liquid and gas flux at the base of the profile were both 0, and heat flux at the boundary was equal to the heat flux sampled by Dakota. Since we were interested in the hydrate saturation near the feather edge, a maximum depth of 1,000 mbsf was used for our model. This ensured that the bottom boundary conditions (constant flux) did not directly influence hydrate and gas predictions at the feather edge.

Dakota and PFLOTRAN Inputs
We used PFLOTRAN to simulate methanogenesis and hydrate formation in a 1D model over time. A period of 120,000 years was chosen, representing the length of a glacial-interglacial cycle. Although this length of time was chosen to represent the glacial-interglacial cycle length, our goal was not to simulate the actual conditions over the last 120,000 years since many of the exact parameters and conditions are unknown to us. Even at the locations we are modeling, there are few constraints on the distribution and amount of gas and hydrate in the system. Therefore, we chose to model a variety of outcomes using Dakota to sample different initial seafloor conditions (TOC, porosity, sedimentation rate, and heat flux) while leaving other conditions such as seafloor temperature and pressure constant over time. This allowed us to get a general sense of how much gas and hydrate could form over a set amount of time which we then used as a basis to predict what might happen if the ocean warms.
To set up this probabilistic model, Dakota was used to provide a distribution of results by sampling the PFLOTRAN input parameters TOC, sedimentation rate, heat flux, and porosity (Table 1). Specifically, Dakota used Latin hypercube sampling to provide a distribution of input variables that were modeled with PFLOTRAN.  Dakota sampled TOC, sedimentation rate, heat flux, and seafloor porosity using a normal distribution. As previously discussed, TOC predictions and standard deviations were modeled using GPSM predictor grids from Phrampus et al. (2020) and observed data points from Seiter et al. (2004) ( Figure 1). To avoid values of zero, a lower bound of TOC was set to 0.01%. An upper bound of TOC was set at 5%, like the GPSM model, to avoid unrealistically high values. Similarly, seafloor porosity was modeled using GPSM predictor grids from Phrampus et al. (2020) and observed data points from Martin et al. (2015). The lower and upper bounds on seafloor porosity were set to 20% and 90% respectively. At depth, porosity was calculated using the trend for marine silty clays presented by Kominz et al. (2011): where ϕ is porosity. Here we assumed the regression of porosity followed the trend presented by Kominz et al. (2011), changing ϕ 0 to the value predicted at a given location in our region of interest. Sedimentation rates and standard deviations were determined from Restreppo et al. (2020) who modeled global oceanic sediment accumulation rates using GPSM at a 5 × 5 arcminute resolution. Bounds on sedimentation rates were set to a minimum of 1 × 10 −14 m/s (3.16 × 10 −5 cm/yr) and to a maximum of 1 m/s. For the studied area, heat flux was sampled from Fuchs et al. (2021), and was set to 48 ± 2 mW/ m 2 . A methanogenesis rate of λ = 5 × 10 −14 s −1 was chosen based on estimates of λ = 1 × 10 −14 s −1 from Bhatnagar et al. (2007) and λ = 1 × 10 −13 s −1 from Malinverno (2010). This value for methanogenesis rate also lies between the constraints of 1 × 10 −15 ≤ λ ≤ 1 × 10 −13 s −1 used by Eymold et al. (2021). We focused our investigation on the area bounded by 36.6°N-38°N and 74.5°W-73.1°W. This area shows some of the higher TOC estimates in the modeled region. Table 2 summarizes the locations where PFLOTRAN and Dakota were used to model hydrate and gas generation. We picked four points along the strike of the continental shelf, located near the feather edge, and five locations along the shelf dip.
At each location, the depth of the seafloor was determined with the Global Multi-Resolution Topography (GMRT) Synthesis (Ryan et al., 2009). The GMRT Synthesis provides high resolution bathymetry data for almost 10% of the global ocean. Individual locations can be queried through the GMRT PointServer Web Service to retrieve accurate water depths at specific latitude and longitude locations.
Seafloor temperatures were calculated at a given depth through a polynomial regression. Data from Boyer et al. (2018) provided 47 temperature profiles from conductivitytemperature-depth (CTD) casts within the 36.6°N-38°N and 74.5°W-73.1°W study area. Since the area of focus was the GHSZ feather edge, near the continental slope, a visual inspection of profile location was used to choose temperature profiles near the continental slope. This left 10 temperature profiles, and a 6 th order polynomial regression was used to create a model of water temperature for depths shallower than 1,000 m: As done with the initial conditions and boundary conditions, the other parameters in the PFLOTRAN model were chosen following the work of Eymold et al. (2021). Many of the values chosen are comparable to values measured at ODP wells in the region. For example, the thermal conductivity parameters for wet and dry sediment of 1 W/m/°C lie between the values found at ODP well site 1,073, ranging between 0.89 and 1.5 W/m/°C (Shipboard Scientific Party, 1998a). ODP wells (sites 902, 903, and 1,073) in this region have permeability measurements in the range of 10 -17 − 10 -16 m 2 (Blum et al., 1996;Dugan et al., 2003). However, we chose to use the permeability value of 10 -15 m 2 as done by Eymold et al. (2021) which can account for siltier beds that may occur in the region.

Geomechanics Model
Following the PFLOTRAN/Dakota simulation of sediment burial and methanogenesis two geomechanical models were used to simulate hydrate dissociation at the BHZ due to an increase in temperature. Our goal was to determine if the conditions needed for sediment failure can exist under the right circumstances, particularly for expected amounts of gas generated microbially, rather than determining the exact conditions needed for sediment failure.
As hydrate dissociates, gas pressure of the system increases causing a reduction in effective stress of the sediment. This can lead to elastic volumetric deformation of the sediment. Although some volumetric expansion will occur due to the increase in gas pressure, shallow marine sediments can still fracture and this fracturing behavior is well described by linear elastic fracture mechanics (Boudreau, 2012;Johnson et al., 2012). Since we are focused near the seafloor, at the feather edge, effective stresses will be small and we are unlikely to build up enough gas pressure to lead to large strains. Therefore, we chose to ignore the volumetric expansion of the sediment. Mechanics for plastic volumetric deformation due to hydrate dissociation have been discussed by Lee et al. (2010) and Waite et al. (2009). In sediments with hydrate saturations of 50-100% Lee et al. (2010) found a decrease in sediment porosity due to hydrate dissociation. However, Waite  (2009) discuss that S h = 25% seems to be a threshold above which pore-filling hydrate becomes load-bearing hydrate. At low saturations, the pore-filling hydrate does not change the shear stiffness of the sediment (Waite et al., 2009). Therefore, in our model we will ignore plastic deformation and possible dilation during dissociation. We focused specifically on the five hydrate bearing locations at the feather edge and modeled the drained and undrained responses that these sediments exhibit during loading. In the drained loading environment, water is allowed to flow out of the sediment during hydrate dissociation and the pore water pressure remains constant. Since the amount of gas formed during dissociation is small, we do not exceed the mobility threshold of gas (Daigle et al., 2020). Therefore, the gas formed will not dissipate through porous flow and will stay in the pore space. In an undrained loading environment, neither water nor gas are allowed to flow and the pore water pressure increases. In both the drained and undrained loading models, the initial pore space of the system was assumed to contain only hydrate and water. The hydrate dissociation reaction was assumed to go to completion (Eq. 4), leading to a final pore space containing only methane gas and water. In addition, no response to dissociation was modeled until the hydrate had fully dissociated to methane gas and water. In this way the final state of the system we model only includes gas and water in the pore space.
In modeling the dissociation reaction only at start and at completion, we ignore changes in pressure within the system that can cause hydrate to reform as local pressure and temperature conditions change during deformation. In addition, we ignore the rate of temperature change and fluid flow as we allow the increase in temperature and final state of the pore space to exist without any additional phase equilibrium calculations.

Drained Model
In the drained loading end member response to hydrate dissociation, water pressure before and after dissociation is constant. Capillary pressure was found using the van Genuchten parameterization: where u c is the capillary pressure, S w is the wetting phase (water) saturation, S wirr is the irreducible wetting phase saturation, P o is the capillary entry pressure, and m is a shape defining parameter (van Genuchten, 1980). The van Genuchten parameters were constrained using mercury intrusion capillary pressure (MICP) measurements performed on marine sediments from locations around the world (Daigle et al., 2020). The following correlations were determined for P o and m: where f c is the clay fraction (mass fraction of particles with equivalent diameter < 3.8 microns) and the reported errors are ±1 standard deviation (Figure 3). Clay fraction was set to 50% (f c = 0.5) to represent a typical marine silty clay and porosity was calculated with depth as described above. To determine S wirr , we used the correlation reported by Daigle et al. (2015): Using the parameters in Eqs 6-8 and f c = 0.5, Figure 3 shows the van Genuchten capillary drainage curve for porosities of 0.5, 0.6, and 0.7. The Peng and Robinson (1976) equation of state was used to calculate the pressure of methane gas based on temperature and the molar volume of the methane gas. We assumed values for methane critical temperature (T c = 190.56 K), methane critical pressure (P c = 45.99 bar), and methane acentric factor (ω = 0.011) based on Poling et al. (2001).
The relationship between gas pressure (u g ) and water pressure (u w ) is given by: Using Eqs 2,5-9 and inputs of water depth, sediment depth, hydrate saturation, original temperature, and final temperature, we solved for capillary pressure using Newton-Raphson optimization. From the capillary pressure solution and the above equations, we were able to calculate final values for water pressure, gas pressure, methane molar volume, gas saturation, water saturation, and capillary pressure.

Undrained Model
In the undrained loading end member response to hydrate dissociation, methane and water mass are both conserved during the dissociation reaction. Thus, mol CH 4 ,i = mol CH 4 ,f and mol H 2 O,i = mol H 2 O,f . To model undrained loading, we used a process similar to the drained loading model. In addition to Eqs 2,5-9, final water saturation (S w,f ) had to be considered as well as water density before and after hydrate dissociation. Water density was determined at a given temperature using the seawater equation of state of Safarov (2003): where u w has units MPa and ρ w is in g/cm 3 . The coefficients A, B, C are polynomial functions of temperature (Safarov, 2003;Safarov et al., 2009). In both the drained and undrained loading responses to hydrate dissociation, porosity was assumed constant. Thus, we were able to use the relation between hydrate and gas saturation from Daigle et al. (2020): where the molar volume of gas is represented by V m,g . We are assuming a structure I hydrate, so the molar volume of hydrate, represented by V m,h , is a constant value: V m,h = 132 cm 3 /mol. In a similar process to the drained loading method, final water pressure was calculated using Newton-Raphson optimization. We then calculated other values using Eqs 2,5-11.

Calculating Stresses
To calculate in situ stresses before hydrate dissociation, a simple calculation for hydrostatic water pressure was used: where z water is water depth in meters, z sed is depth below seafloor in meters, ρ w = 1,024 kg/m 3 , g = 9.81 m/s 2 , and u w is in Pa.
Overburden stress (σ v ) was determined by: where σ v is in Pa. Bulk density (ρ b ) was determined at a given depth using the calculated porosity curve: with a grain density ρ g = 2,700 kg/m 3 . Effective stress (σ′) prior to hydrate dissociation can be simply calculated as: Following the methods of Daigle et al. (2020), we assumed that the maximum principal stress is vertical (σ 1 = σ v ) and the sediments are vertically transversely isotropic (σ 2 = σ 3 = σ h ). The relation of σ h ′ and σ v ′ can be expressed as: or where ] is Poisson's ratio and subscripts v and h represent vertical and horizontal stresses respectively. To calculate Poisson's ratio with depth, a sixth-order polynomial was fit to marine sediment data up to 650 m deep from Hamilton (1979).
After dissociation, the pore space is a two-phase system of methane gas and water. To calculate effective stress, we followed the methods of Bishop (1959) and Nuth and Laloui (2008): where the equivalent pore pressure, π, is a weighted average of the water pressure and gas pressure in the system. Eqs 19,20 also hold for the pre-dissociation case with no gas pressure as u g = 0 and thus π = u w . With this, the change in effective stress at the BHZ due to hydrate dissociation was calculated.

Failure Model
A nonlinear Hoek-Brown failure envelope was constructed at each location with two goals in mind: 1. to determine if any failure occurred due to hydrate dissociation, and 2. which type of failure occurred (tensile or shear failure). This followed the work of Hoek and Brown (1997) and the methodology of Daigle et al. (2020). We assumed that the sediments were initially intact with no faults, fractures, or joints, so the tensile strength, T, was calculated by: where c u is the unconfined compressive strength and m i is the Hoek-Brown constant. We set m i = 4 as reported by Hoek and Brown (1997) for claystones. This also follows the work done by To calculate c u , we used the correlation derived for muds and mudrocks by Ingram and Urai (1999): log 10 c u −6.36 + 2.45 log 10 0.86V p − 1172 (22) where c u is in MPa and V p is the compressional wave velocity in m/s. At each depth we calculated V p using the equation derived from Hamilton (1979) for marine sediments: Finally, the Hoek-Brown failure envelope was calculated from Hoek and Brown (1997): where τ is the shear stress of the sample, σ n ′ is the effective normal stress of the sample, and A and B are Hoek-Brown fitting parameters. Values of A = 1.02 and B = 0.858 were calculated by fitting failure envelopes to the results of triaxial tests performed on marine clays from the Gulf of Mexico and the Bay of Bengal (Silva et al., 2000;Moses et al., 2003;Dugan and Germaine, 2009). The Hoek-Brown failure envelope was calculated based on the initial sediment properties prior to dissociation. Once dissociation occurred, both σ v ′ and σ h ′ were recalculated using Equation 19. Failure occurs if any points of the new Mohr circle lie to the left of the Hoek-Brown failure envelope. Since we assumed that the dissociation reaction went to completion, it was possible that the new effective stresses (σ v,f ′ and σ h,f ′ ) both lie past the Hoek-Brown failure envelope. In reality, failure would initiate as soon as the Mohr circle began to cross the Hoek-Brown envelope: tensile failure occurring if the Mohr circle crossed the failure envelope at σ 3 ′ T (where σ 3 ′ σ h ′) and shear failure occurring if the Mohr circle crossed the failure envelope when |σ 3 ′| < |T| (Figure 4).

GPSM
The weight percent of TOC at the seafloor predicted with GPSM is shown in Figure 5 as well as the standard deviation and mean inexperience of these predictions. The highest TOC values were predicted along the line from (35.4°N, 75.0°W) to (39.0°N, 72.0°W). Comparing this region in the prediction map to the standard deviation map, known seafloor TOC values are sparse in this area and TOC standard deviations are higher than the rest of the predicted grid. The locations chosen along the shelf strike can be found in this area. The mean inexperience is calculated as the average distance in parameter space from the predicted location to its five nearest neighbors. Compared to the standard deviation map, a relatively low mean experience can be seen along the line from (35.4°N, 75.0°W) to (39.0°N, 72.0°W). Thus, although there is a high variance of predictions, there are other locations globally that are parametrically close to locations along this shelf strike line. A summary of the predicted weight percent of TOC for the nine modeled locations can be found in Table 3.

Dakota and PFLOTRAN
The range of input values sampled by Dakota and used in the PFLOTRAN burial process to generate gas and hydrate are summarized in Table 3.
Average TOC predictions were around 2% dry weight at the shallower locations near the feather edge (sites 1-6). At deeper locations (sites 7-9), TOC predictions were closer to 1% dry weight. Seafloor porosity was around 65% at shallow locations and predictions increased to 73% for deeper location. Average sedimentation rates were just over 0.007 cm/yr with standard deviations around 0.1 cm/yr (Restreppo et al., 2020). As discussed in the Methods section, to avoid negative sedimentation rates, a lower bound for sedimentation rate was set at 3.16 × 10 −5 cm/yr and any values sampled smaller than this value were set to 3.16 × 10 −5 cm/yr. These sedimentation rates are lower than the sedimentation rates on the U.S. Atlantic margin found to the north and south, but the studied location could simply be in an area of low sediment influx. Measurements to the north along the U.S. Atlantic margin (offshore of New Jersey) range from 0.011 cm/yr to 0.058 cm/yr (Shipboard Scientific Party, 1994a;Shipboard Scientific Party, 1994b;Shipboard Scientific Party, 1994c;Shipboard Scientific Party, 1994d;Shipboard Scientific Party, 1994e;Shipboard Scientific Party, 1998a) while locations to the south (offshore North Carolina) are closer to 0.02 cm/yr (Shipboard Scientific Party, 1994b;Shipboard Scientific Party, 1994c;Shipboard Scientific Party, 1996a;Shipboard Scientific Party, 1998b). Location nine had a much lower predicted sedimentation rate than the other eight locations, corresponding to its further distance from the shelf edge and seacoast, as well as its overall depth. Heat flux for all simulations were between 42.3 and 52.2 mW/m 2 . At each location, we ran 50 simulations using PFLOTRAN, and sampled the seafloor TOC, seafloor porosity, sedimentation rate, and heat flux inputs with Dakota. For an individual location, the resulting output profiles of gas, hydrate, and temperature were plotted against depth, and a base of hydrate stability was calculated. Figure 6 shows the hydrate saturation profiles and the gas saturation profiles with depth for the nine study locations. Figure 7 shows an example temperature profile with depth for site 1, and includes the hydrate phase envelope calculated from Kamath (1984): where hydrostatic pressure, P (kPa), can be calculated at a given sediment depth using Eq. 12, allowing us to calculate temperature, T (K). The phase envelope was included to illustrate the difference in the BHZ between simulations. Due to the diversity of sampled heat flux values, the geothermal gradient for each simulation is different, leading to a range of depths where the phase envelope is crossed and thus a range of BHZ depths. Table 4 summarizes the hydrate profiles shown in Figure 6. Some simulations failed to run or ran into oscillatory errors, thus producing no hydrate or gas predictions. The number of runs that produced hydrate and gas profiles (successful runs) is noted at each location. In addition, the number of runs where hydrate formed is summarized. At each site, the average BHZ was calculated as the depth where the hydrate phase envelope from Kamath (1984) intersected the average temperature profile at the location. In addition, we noted the maximum depth at which hydrate formed (Max hyd. depth) at each location. This depth represents the extent of the zone of actual hydrate occurrence where the concentration of methane exceeds the solubility of methane in seawater in addition to the pressure and temperature conditions of the hydrate stability zone (Xu and Ruppel, 1999). At sites 1-4 and site 6, the average extent of the zone of hydrate occurrence is very similar the calculated BHZ. Looking at the profiles for these five locations, hydrate concentrations are often at a maximum at this depth before they drop to zero and methane gas is witnessed instead. The minimum and maximum extent of the zone of hydrate occurrence are also noted, illustrating the varying extent of the BHZ between models as shown in Figure 7.
At sites 7-9, fewer simulations result in the growth of hydrate. As was done for sites 1-4 and 6, the average depth of the BHZ was calculated for these sites. Interestingly, at these deeper locations, the average extent of the zone of hydrate occurrence was much shallower than the depth of the BHZ. At site 7, the average BHZ is 516.5 mbsf, although graphically it can be seen that a majority of the hydrate profiles do not reach this depth ( Figure 6). Even fewer simulations form hydrate at sites 8 and 9, and none of the simulations produce hydrate near the depth of the BHZ (601.4 mbsf and 647.7 mbsf respectively) at these locations. In profile 5, no hydrate forms. This location is upslope of the feather edge and outside of the hydrate stability zone.
A distribution of total hydrate mass for each profile is shown in Figure 8. Site 6 has the largest distribution of total hydrate formed, however site 7 has the largest amounts of hydrate formed in a run. The larger values of hydrate mass in sites 7-9 can be attributed to the thickness of the hydrate stability zone at these locations. Near the feather edge (sites 1-4) many runs do not generate hydrate, and those that do have very thin hydrate layers due to the shallow BHZ. This leads to much lower predictions of total hydrate formed in this area.
The gas saturation profiles shown in Figure 6 are summarized in Table 5. Gas generation occurs within the modeled 1,000 m sediment column in eight of the simulations. The only simulation in which no gas was generated during the 50 runs was at the deepest location, site 9. When calculating the maximum gas saturation and the depth of the maximum gas saturation in Table 5, only maximum saturations shallower than 1,000 m were considered. Therefore, site eight does not have a maximum gas saturation value or depth. The maximum gas saturation for site seven occurs in only one run and may be an outlier. Comparing the remaining six locations, site 5 has the largest average maximum gas saturation. As discussed, this site is upslope of the feather edge and outside of the hydrate stability zone.

Dissociation and Failure Criteria
At the five locations along the feather edge (sites 1-4, 6) we modeled hydrate dissociation at the BHZ due to a temperature increase of 1°C, ignoring any changes in pressure over time due to the rise and fall of sea level were ignored. At each location, we calculated a final gas saturation due to dissociation, along with the resulting final water pressure, final gas pressure, and capillary pressure ( Table 6). A value for π (Eq. 20) was also calculated at Frontiers in Earth Science | www.frontiersin.org May 2022 | Volume 10 | Article 835685 each location. Comparing the drained and undrained models, both predict a decrease from the initial saturation of hydrate to the final saturation of gas. This is because a constant density (and thus molar volume) was chosen for hydrate whereas the methane gas that forms is compressible. Although the BHZ at these locations is relatively near the seafloor, there is still enough pressure from the 50-100 m of sediment and 500 m of water to compress any gas that forms. In the drained model, hydrate saturation decreases to the final gas saturation by a factor of about four ( 1 4 S h ≈ S g ). In the undrained model, the decrease in hydrate saturation to gas saturation is over twice as large ( 1 10 S h ≈ S g ). In the undrained environment water pressure is not held constant during dissociation and instead the total moles of water before and after dissociation is constant, leading to a larger decrease in hydrate to gas saturation. The impact of drainage can also be seen when comparing the water pressures calculated before (u w,i ) and after (u w,f ) dissociation for both models. By definition, the drained model has a constant water pressure during hydrate dissociation. This leads to final gas pressures (u g,f ) that are similar to the water pressure. On the other hand, when hydrate dissociates in the undrained model, water pressure nearly doubles at every location. Capillary pressure (u c ) in the undrained model is similar to capillary pressure in the drained model, however the large increase in water pressure leads to a higher final gas pressure and higher π value (Eq. 20).
At each location, the total horizontal and vertical stresses and the effective horizontal and vertical stresses were calculated before and after hydrate dissociation ( Table 6). Total stresses were relatively consistent among the five locations as they all lie along the shelf strike. During dissociation, the effective stress on Hydrates are stable where the temperature profile is to the left of the hydrate phase envelope (black dashed line). Around 80 mbsf sediment depth at Site 1 (water depth 508 m), temperature profiles begin to cross the hydrate phase envelope, so hydrate is no longer stable. At the deeper Site 8 (water depth 2,595 m), temperature profiles begin to cross the hydrate phase envelope around 550 mbsf. Note the different depth scales between the two temperature profiles. This was done to better visualize the intersection of the phase envelope at Site 1. n/a n/a n/a n/a n/a n/a   n/a n/a 9 50 50 0 n/a n/a

GPSM
Our predictive map of seafloor TOC from GPSM was comparable to the work done by Lee et al. (2019) and Eymold et al. (2021). The global map of seafloor TOC produced by Lee et al. (2019) also predicted an area of increased TOC off the east coast of the United States in the same region [(35.4°N, 75.0°W) to (39.0°N, 72.0°W)] predicted in our model. The map of seafloor TOC predicted by Eymold et al. (2021) covers about a quarter of the area of our prediction. In this prediction, trends of TOC are similar to the trends we predict as Eymold et al. (2021) also predicts an increase in TOC at the southern end of the (35.4°N, 75.0°W) to (39.0°N, 72.0°W) region. It is unsurprising that the areas with increased TOC predictions are similar between this work and the work done by Lee et al. (2019) and Eymold et al. (2021). Although the data sets used between the three predictions differed slightly, the general methodology was consistent across the models.
To ensure the accuracy of our seafloor TOC model, it is additionally important to compare the results to field data in the area. Investigations by Skarke et al. (2014) have found multiple instances of methane gas leakage from the seafloor along the U.S. Atlantic margin in the form of gas plumes (Figure 9). The sites of these methane seeps correspond with areas predicted by our GPSM model to have higher values of seafloor TOC, specifically in the region between (35.4°N, 75.0°W) and (39.0°N, 72.0°W). The seeps identified by Skarke et al. (2014) are mostly found just updip of the feather edge. At these locations, seepage is most likely due to small tensile failures that occur when microbial gas is generated close to the seafloor (Daigle et al., 2020). However, at the feather edge where we modeled hydrate and gas formation with PFLOTRAN/Dakota, if hydrate dissociates it will instead cause shear failure.
Around the U.S. Atlantic margin, the U.S. Bureau of Ocean Energy Management (BOEM) identified six unique regions with BSRs of moderate and high confidence (Frye et al., 2013). Some of these regions with identified BSRs overlay the region where we predicted high TOC values and modeled hydrate and gas formation. Our model forms hydrate and gas through methanogenesis and therefore suggests that the BSRs identified by Frye et al. (2013) may have a similar local microbial origin and that the concentration of gas and hydrate in the pore space may be small.

Hydrate and Gas Formation
When using PFLOTRAN, a 120,000 years time period corresponding to the length of a glacial-interglacial cycle was chosen to simulate methanogenesis and hydrate formation. This gave us a good idea of hydrate profiles along the feather edge of hydrate stability, and we demonstrated that conditions exist in which failure will occur due to dissociation. If we instead wanted to model the actual burial history and hydrate formation in this area, we would need to include changes in temperature, pressure, and sedimentation rate over time. However, even though gas seeps (Skarke et al., 2014) and regions with BSRs of moderate to high confidence (Frye et al., 2013) have been identified in the area we are modeling (Figure 9), there are still few constraints on the distribution and amount of gas and hydrate in the area. Therefore, we chose to model a variety of outcomes using Dakota to sample different initial seafloor conditions.
Once the hydrate/gas profile has been put in place we assumed an increase in temperature on a time scale much shorter than 120,000 years. One cause of this temperature change could be a weakened Atlantic Meridional Overturning Circulation (AMOC) which could lead to increased surface and seafloor temperatures in our area of interest (Liu et al., 2020;Garcia-Soto et al., 2021). The movement of the feather edge due to changes in seafloor temperature has also been shown on a shorter timescale by Phrampus and Hornbach (2012) just south of the region we modeled. Changes in sea level can also alter the GHSZ as pressure changes as the seafloor. Not only can changes in sea level shift the feather edge, but lower past sea levels can create an area within the FIGURE 9 | Output maps from GPSM over the area 32°N-41.5°N and 78°W-68.5°W. Methane seeps (yellow dots) identified by Skarke et al. (2014) are also plotted and support the increased TOC prediction between (35.4°N, 75.0°W) and (39.0°N, 72.0°W). Locations where we modeled hydrate and gas formation are marked with a magenta × and numbered. Red and yellow outlines are areas where the U.S. Bureau of Ocean Energy Management (BOEM) identified BSRs in seismic data (Frye et al., 2013). current GHSZ where hydrate was unable to accumulate over time (Stranne et al., 2016).
In the area of study, varying sedimentation rates over the Pleistocene are noted by McHugh and Olson (2002). Since our model relies on sedimentation to transport organic matter, it would seem that given identical TOC measurements, the sedimentation rate would increase the amount of methane hydrate and gas formed in the simulation. This is discussed by Eymold et al. (2021) who concur that as the sedimentation rate approaches zero, hydrate growth is unexpected. Eymold et al. (2021) also note that increasing sedimentation rate in their experience can eventually lead to reduced hydrate growth as the organic carbon is buried so quickly it cannot go through the methanogenesis process. In reality, changing sedimentation rates over time will affect the concentration of TOC on the ocean floor.
In the area we model, organic matter on the seafloor is mainly marine-driven with only a little mixing of terrestrial material (Shipboard Scientific Party, 1994a;Shipboard Scientific Party, 1994b;Shipboard Scientific Party, 1994c;Shipboard Scientific Party, 1994d;Shipboard Scientific Party, 1994e). Higher sedimentation rates in the area would lead to a larger flux of terrestrial organic to the seafloor, however there would also be more sediment deposited relative to the standard background flux of marine organic carbon. Therefore, a variable sedimentation rate over time would have to consider the changing seafloor TOC concentrations and may still result in predictions similar to models where seafloor TOC and sedimentation rates are held constant.
Sampling methanogenesis rate in the PFLOTRAN/Dakota model also changes the outcomes of the simulations. Since we do not know the exact methanogenesis rate in the area or the probability distribution of the methanogenesis rate, we chose a constant value within the range modeled by Eymold et al. (2021) and also between the values of λ estimated by Bhatnagar et al. (2007) and Malinverno (2010). The effect of changing methanogenesis rate in the model is discussed by Eymold et al. (2021) who found that higher methanogenesis rates led to a deeper BHZ but a lower maximum hydrate saturation. We have included additional hydrate and gas profiles in the supplementary materials where we sampled methanogenesis rate in addition to the other previously sampled parameters.

Sediment Failure Model
In all five locations along the feather edge (sites 1-4, 6) shear failure occurred in the undrained model while no shear failure occurred in the drained model. At these sites, the largest average saturation of hydrate at the BHZ was predicted to be 1% (site 6). At the other four sites, the average predicted hydrate saturation at the BHZ was at or below 0.7%. At these locations gas saturation after dissociation was no more than 0.2%. In the undrained environment, gas saturations formed from dissociation were below 0.1%. Even with these relatively low saturations of gas formed from dissociation, hydrate dissociation under undrained conditions led to shear failure. This differs from the conclusion of Daigle et al. (2020) who suggest that tensile fracturing is generally favored in the shallowest sediments near the feather edge. The BHZ modeled in this paper with PFLOTRAN is at depths between 40 and 80 m (along with site 6 at 182 m). Thus, the sediments modeled may still be too deep below the seafloor to experience gas-driven tensile fracturing.
In the undrained model, shear failure occurs due to increased pore pressure and gas pressure in the sediment. As pore pressure increases, the effective stress decreases to the point where the sediment enters a failure regime. As depicted in Figure 4, this can be visualized as the Mohr circle shifting to the left until it intersects with the Hoek-Brown failure envelope. Due to the depths of our modeled BHZ, the representative Mohr circle intersected with Hoek-Brown failure envelope at |σ 3 ′| < |T|. Closer to the seafloor the stresses are nearly isotropic, and the representative Mohr circle could be small enough that a decrease in effective stress leads to tensile failure when |σ 3 ′| |T|. Stranne et al. (2017) note that fractures will not occur near the feather edge in sediments with permeabilities larger than 10 -15 m 2 . Although we use a permeability value of 10 -15 m 2 in our burial model to account for siltier beds, multiple ODP wells (sites 902, 903, and 1,073) in the region describe permeabilities between 10 -17 and 10 -16 m 2 (Blum et al., 1996;Dugan et al., 2003), so shear failure still may occur due to hydrate dissociation.
The low saturation of gas formed from hydrate dissociation needed for the dissociation to cause shear failure has interesting implications on slope stability for the region. Looking at the area of high TOC values predicted with GPSM that overlap with the methane seepage locations highlighted by Skarke et al. (2014), there is a 250 km region along the continental slope where high TOC values are expected (and where methane seeps have already been found). Even with low hydrate concentrations between 0.2% and 1%, if this hydrate dissociates in an undrained environment, shear failure could ensue. In reality, the drainage environment will be somewhere between the drained and undrained endmembers. Therefore, it will eventually be important to have a mixed-drainage model to determine the possibility of failure due to hydrate dissociation in a region. In addition, it will be important to look at hydrate dissociation as a step by step process where local pressure change during dissociation can cause hydrate to reform. When considering this, the rate at which temperature changes and the rate of fluid flow are important in correctly modeling the grain space.
For slope failure or submarine landslides to occur in an area, there needs to be some sort of triggering mechanism such as an earthquake, changing sea level, or high internal pore pressure causing the downslope component of stress to become larger than the resisting stress (Hampton et al., 1996). The dissociation of methane hydrate has also been found to lead to shear failure in seafloor sediments resulting in seafloor slumping and submarine landslides (Kvenvolden, 1993;Hampton et al., 1996;Hornbach et al., 2007;Elger et al., 2018). The timing of slope failure due to hydrate dissociation has been investigated by Lee (2009), Maslin et al. (2004, and Nixon and Grozic (2007), who have shown possible connections globally between hydrate dissociation and submarine landslides. The mechanism by which hydrate dissociation leads to slope instability has been discussed as well. In both fine grained and coarse grained sediments, fluid migration towards the seafloor can occur along active thrust faults leading to a build up of gas pressure, resulting in lower effective stresses and sediment instability (Conti et al., 2008;Argentino et al., 2019). In areas without thrust faults to act as a conduit for gas migration, Kvenvolden (1993) suggests that the dissociation of hydrates can create an "enhanced fluidized layer" beneath the GHSZ, leading to the slope failure. Elger et al. (2018) suggest that the overpressure of gas at the BHZ creates a gas pipe to shallower sediments, and it is in these shallower coarse-grained sediments where gas builds up, leading to shear failure.
A factor of safety calculation, containing the seafloor slope variable, compares resisting to shearing forces in a sediment and can be used to determine if slope failure is likely to occur at a submarine location (Løseth, 1999;ten Brink et al., 2009). Along the U.S. Atlantic continental slope most submarine landslides originate from areas with a seafloor slope angle between 2 and 4° (Booth et al., 1993). In the area of our study, the gradient of the continental slope is mostly between 4 and 8 with some areas reaching 8-12° ) thus slope failure may be a concern under the undrained conditions that we modeled.

CONCLUSION
We used a geospatial machine learning model to create a global map of TOC on the seafloor. Focusing specifically on the U.S. Atlantic margin, the region with high TOC predictions between (35.4°N, 75.0°W) and (39.0°N, 72.0°W) was consistent with methane gas seeps located by Skarke et al. (2014). We then modeled hydrate and gas formation over a 120,000 years time period at nine locations in this area: five hydrate bearing locations along the feather edge, one gas location updip of the feather edge, and three hydrate bearing locations downdip of the feather edge. For each of these locations, hydrate and gas formation were modeled by sampling seafloor TOC, seafloor porosity, sedimentation rate, and heat flux.
At the feather edge, we modeled hydrate dissociation at the BHZ due to an increase in temperature while ignoring the change in pressure due to the rise and fall of sea level. In a purely drained environment, no failure is expected to occur due to hydrate dissociation. However, in an undrained environment, the criterion for shear failure is quickly met during hydrate dissociation. Even as gas saturations due to hydrate dissociation stayed below 0.1% (hydrate saturations: 0.2%-1%) shear failure was predicted to occur. This suggests that the hydrate that forms over one glacial cycle can cause submarine slope failure upon dissociation.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.