A New Dynamic Modeling Approach to Predict Microbial Methane Generation and Consumption in Marine Sediments

: Methane, as a clean energy source and a potent greenhouse gas, is produced in marine sediments by microbes via complex biogeochemical processes associated with the mineralization of organic matter. Quantitative modeling of biogeochemical processes is a crucial way to advance the understanding of the global carbon cycle and the past, present, and future of climate change. Here, we present a new approach of dynamic transport-reaction model combined with sediment deposition. Compared to other studies, since the model does not need the methane concentration in the bottom of sediments and predicts that value, it provides us with a robust carbon budget estimation tool in the sediment. We applied the model to the Blake Ridge region (Ocean Drilling Program, Leg 164, site 997). Based on seaﬂoor data as input, our model remarkably reproduces measured values of total organic carbon, dissolved inorganic carbon, sulfate, calcium, and magnesium concentration in pore waters and the in situ methane presented in three phases: dissolved in pore water, trapped in gas hydrate, and as free gas. Kinetically, we examined the coexistence of free gas and hydrate, and demonstrated how it might affect methane gas migration in marine sediment within the gas hydrate stability zone.


Introduction
Marine sediments are the largest known methane reservoir on Earth [1][2][3][4]. The shortterm carbon cycle is linked to the long-term geological sequestration of carbon by the organic matter mineralization/degradation and preservation in marine sediments [5,6]. One of the best ways to understand the global carbon cycle and its effects on, for instance, climate change, is to quantify the rates of biogeochemical processes in marine sediments [7]. The rate of organic matter degradation decreases with sediment depth, age [8], and quantity of metabolites [9].
A significant mass of CH 4 is stored in marine gas hydrates, which are a potential energy resource [10] but also able to impact the atmospheric inputs of methane and, therefore, affect Earth's climate [11][12][13]. Gas hydrate also affects the global carbon cycle since they hold vast quantities of carbon (~550 Gt of carbon) [14,15] and trap large amounts of free gas below the gas hydrate stability zone (GHSZ) [16]. As methane in the gas hydrate is dominantly microbial in origin [17], the study of organic matter decomposition kinetic in anoxic sediments has increased [1,18,19]. Numerical modeling can improve the insights into the physical and chemical processes of hydrate formation, which is a massive sink of carbon in sediments [20,21]. According to the thermodynamic perspective, free gas should not coexist with pore water within the GHSZ except along the triple point line [22][23][24]. Nevertheless, widespread seafloor methane venting reported in many regions of the world suggests the coexistence of free gas and hydrate in those places [25][26][27][28][29][30][31][32][33][34]. Thermal anomalies [35], salinity [36], capillarity [37], and crustal fingering [38] are some explanations for these field observations. Although many studies have addressed this subject from the thermodynamic point of view, few efforts have been devoted to understanding gas hydrate systems kinetics, which facilitates the free gas flow through GHSZ.
Since marine sediments are the major component of the global carbon cycle [2,4] and they have significant impacts on the oceans and sediment chemistry, and Earth's climate, this study has focused on dynamic modeling of sedimentary organic matter mineralization and the generation and consumption (via sulfate reduction) of methane in sediments. Data obtained from the Ocean Drilling Program (ODP), in the Blake Ridge region, in the North Atlantic U.S. margin (Leg 164, site 997) was used to validate the model. The current study's approach is to use a new numerical model to enhance the general understanding of biogeochemical processes by considering progressive sedimentation over time and provide a more reliable dynamic approximation of processes in the marine sediments.
In contrast with other studies [1,4,6,9,12,19,[39][40][41][42][43][44], which used the data of seafloor and bottom of sediments, the major advantages of our transport-reaction model are using the seafloor data as the only input and simulating sedimentation over time with progressive compaction of sediments. Since most of the time, it is hard or impossible to obtain the data in the deeper part of sediments, in this model, we provide a unique tool to predict methane in all phases (pore water, gas hydrate, and free gas), in addition to the concentration of sulfate, dissolved inorganic carbon, calcium, and magnesium in pore water. Furthermore, applying our model in different areas makes it possible to have a reliable prediction of the amount of methane formed in sedimentary basins, even with limited available data. Moreover, by changing the input data (such as POC and sedimentation rate) with time, it is possible to investigate their effect on biogeochemical processes. Additionally, here, we address the free gas and hydrate coexistence within the GHSZ from the kinetic perspective, which could be a new direction for future research in this area.
The rest of the paper is organized as follows. In Section 2, data source and numerical modeling of solid, dissolved, and gas species are explained. The results of simulation and comparison with available data as well as discussion are presented in Section 3. The final section summarizes the main conclusions obtained through this study.

Data Source
Modeling was done using available data from Janus database (http://www-odp.tamu. edu/database/), which stores data for ODP. According to the excellent pore water data set, we applied our model on site 997, Leg 164 (Blake Ridge, North Atlantic U.S. margin).
Site 997 was drilled at the Blake Ridge crest, a large contourite system [45] located in the North Atlantic U.S. margin ( Figure 1). This site was selected to validate our numerical model results due to the availability of high-resolution pore water chemistry data, the presence of gas hydrate, and abundant methane of microbial origin [46]. Well-developed bottom simulating reflections (BSR) are recognized in seismic sections in the area and are interpreted as the result of free gas accumulation below the GHSZ [47]. The drilled deposits consist of Holocene to late Miocene, light gray and greenish gray, fine-grained (clay) sediment with abundant (up to 50 wt%) authigenic carbonate, which occurs as nodules and concentrated beds mostly between 50-100 m below seafloor (mbsf) and >500 mbsf (Figure 1). The occurrence of dispersed gas hydrate was reported between 180-450 mbsf, and a large piece (15 cm long) was recovered at 331 mbsf.  [48]). (b) Composed core-log for site 997 showing calcium carbonate (wt%) in sediments. Ellipses and bars indicate the position of carbonate nodules and carbonate beds, respectively (modified from Paull et al. [49]). The blue area shows the depth range for the occurrence of dispersed gas hydrate in sediments, and the star shows the location of a 15 cm long gas hydrate piece.

Numerical Modeling
A dynamic numerical transport-reaction model [50] is applied to simulate particulate organic matter mineralization in marine sediments over time, under anoxic conditions. In this study, sediment deposition with its physical properties is reconstructed over time using sedimentation rate, geothermal gradient (Table 1), and compaction. The model calculates the spatial and temporal evolution of ten variables, i.e., concentration versus depth profiles of four solid species: (1) particulate organic carbon (POC), (2) calcium carbonate (CaCO 3 ), (3) magnesium carbonate (MgCO 3 ), and (4) gas hydrate (GH); five dissolved species: (5) sulfate (SO 4 2− ), (6) methane (CH 4 ), (7) dissolved inorganic carbon (DIC), (8) calcium (Ca 2+ ), and (9) magnesium (Mg 2+ ), and one gas species: (10) methane. Organic matter degradation, sulfate reduction, methanogenesis, anaerobic oxidation of methane (AOM), carbonate precipitation, and methane hydrate formation and dissociation are the main processes considered in the model. The model solves the partial differential equation for dissolved species (1), solids (2), and gas (3) based on the classical approach used in early diagenesis modeling [7,9,52] where t is time, x is depth in the sediment, C l is the concentration of dissolved species in pore water, ∅ is the porosity, D s is the molecular diffusion coefficient in pore water corrected for tortuosity, sediment temperature, salinity, and pore water viscosity, v is the burial velocity of solutes, C s is the concentration of solids in dry sediments, w is the burial velocity of solids, C g is the concentration of free gas (normalized to pore water), D g is the apparent gas diffusion coefficient, v g is the relative velocity of ascending gas, and R stands for the reactions occurring in the sediment column. The decrease in porosity with sediment depth, advective transport of solids and solutes via burial, and steady-state compaction are considered in the model. The model assumes negligible bioturbation and bioirrigation since these are typically confined to the upper 10-20 cm [1].
Our model consists of a column (1D) extending from the seafloor down to a depth of 600 mbsf. Constant concentrations of dissolved species and solids were considered as the upper boundary of the column (Dirichlet boundary conditions), corresponding to the seafloor, based on their seawater concentrations. According to ODP data for POC in site 997, a function was defined for POC input with time. POC input is considered 1.6 wt% and decreases to 0.56 wt% during the last 1.5 Myr. Zero concentration gradients (Neumann boundary condition) were prescribed as the lower boundary condition (base of the model at 600 mbsf) for sulfate, calcium, and magnesium. According to the literature, there are significant gradients of methane and DIC at lower boundary conditions, indicating that dissolved species are released in deeper sediments and transported to the surface by diffusion [9]. For simplification, in most cases, they considered a constant value or zero concentration gradients at the bottom of their modeled columns [1,9,15,43,44,53]. In this study, according to the depth and time, two boundary conditions were considered because of the sedimentation evolution. According to time, before the first layer of the sediment reaches the lowest point of the modeled column, a zero concentration gradient is imposed for methane, DIC, and free gas. After that time, by extrapolating the last two points of DIC, methane, and free gas, a virtual point was considered below the last point and used for the lower boundary condition in the bottom of the modeled column ( Figure 2). Therefore, these values were changing over time to have a more realistic and novel dynamic model. A detailed description of the model and the list of used parameters are given in Appendix A.
Simulating the growing sediment column over time with progressive compaction of sediments is a significant advantage of this model, so the model was run for 100 million years to simulate an extended period.
The ten sets of the differential equations consist of four solid species (using Equation (1)), five dissolved species (using Equation (2)), and one gas species (using Equation (3)), and were discretized using finite difference techniques and the method-of-lines [54]. The model was solved using MATLAB (R2020a, MathWorks). All modeling and calculations were done based on the data gathered from OPD site 997, Leg 164 [22,49,51]. sediment reaches the lowest point of the modeled column, a zero concentration gradient is imposed for methane, DIC, and free gas. After that time, by extrapolating the last two points of DIC, methane, and free gas, a virtual point was considered below the last point and used for the lower boundary condition in the bottom of the modeled column ( Figure  2). Therefore, these values were changing over time to have a more realistic and novel dynamic model. A detailed description of the model and the list of used parameters are given in Appendix A.

Solid Phase and Pore Water Concentration
The results in the last time step of the model are shown in Figure 3. The timelapse video of the model evolution is presented in Video S1. According to the results, sulfate concentration decreases linearly from the seafloor to the depth of 24 mbsf, which is the depth of the sulfate-methane transition (SMT). Between 0 and 18 mbsf, methane concentration is almost zero; then, its concentration rapidly increases, indicating ongoing methanogenesis and the termination of methane consumption via anaerobic oxidation of methane (AOM) in the sulfate reduction zone. In the methane profile (Figure 3b), the quantity of methane generated reaches saturation below about 170 mbsf. The higher methane values measured in the ODP reference site 997 suggest a possible dissociation of gas hydrate and methane release during core recovery.
Considering the seawater as the only source of calcium and magnesium in the system, the model predicts a decline in magnesium and calcium concentrations with depth due to carbonate precipitation associated with organic matter mineralization and related alkalinity increase at shallow depths below seafloor (<100 mbsf) [55]. Despite a slight increase in calcium content at about 500 mbsf in the study site, the model shows a good match with the reference site's data. An additional source of calcium in the sediment and change in abundance of authigenic carbonate [51] could explain the difference between the model and data.
In the model, DIC results from the production of sulfate reduction, AOM, and methanogenesis [9]; however, DIC data presented here for ODP is the sum of total alkalinity, mostly produced by sulfate reduction and AOM [56], and carbon dioxide produced by methanogenesis. depth of the sulfate-methane transition (SMT). Between 0 and 18 mbsf, methane concentration is almost zero; then, its concentration rapidly increases, indicating ongoing methanogenesis and the termination of methane consumption via anaerobic oxidation of methane (AOM) in the sulfate reduction zone. In the methane profile (Figure 3b), the quantity of methane generated reaches saturation below about 170 mbsf. The higher methane values measured in the ODP reference site 997 suggest a possible dissociation of gas hydrate and methane release during core recovery.  total alkalinity + CO 2 ) in pore water, and (d) weight percent of particulate organic carbon (POC) in sediment. Lower row: concentration of (e) calcium (Ca 2+ ); (f) magnesium (Mg 2+ ) in pore water, and the volume percent (calculated based on the pore volume) of (g) hydrate (GH), and (h) free gas (CH 4 (G)). Figure 4 shows methane generation and sulfate consumption (via sulfate reduction and anaerobic oxidation of methane) simulated over time. At the top of the model (0-23 mbsf), most organic matter is consumed by the sulfate reduction, and the concentration of methane is nearly zero at all times. Below that depth, methane concentration slightly increases because of the onset of methanogenesis and cessation of sulfate reduction owing to sulfate depletion. After about 7 million years, more methane is produced, and the rate of AOM increases to reach an almost constant and relatively steady-state situation. As shown in Figure 4, at the beginning of the model (0-0.34 Myr), the SMT occurs in a deeper position (42 mbsf) in relation to present day (24 mbsf) because there is no methane below to be consumed via AOM which would compete for sulfate with sulfate reduction fueled by POC. As time goes on and more methane is produced, the diffusion of methane to the upper part increases and results in faster depletion in sulfate by AOM, with consequent shallowing of the SMT depth (20 mbsf). According to the decrease in POC input (98.5-100 Myr) based on the reference data for site 997, methanogenesis and sulfate reduction slightly decrease, and SMT depth increases to its present day (24 mbsf).  According to the ODP report for the studied site [51], Blake Ridge contains gas hydrate, and below the GHSZ, there is an accumulation of free gas [16,57], which is produced via microbial methanogenesis [22]. According to the salinity profile, we considered up to 8% of pore volume could be occupied by gas hydrate and free gas. As shown in Figure 4, based on the model, methane produced by methanogenesis are in three phases, dissolve in pore water, as gas hydrate, and free gas. In the beginning, there is no methane in the system. By consuming all sulfate in depth, methane starts to be produced by methanogenesis. Over time, more methane is generated, and the concentration of methane in pore water increases, so it diffuses to the upper part of the column. After around three million years, when about 350 m of sediment accumulated, methane in pores below about 180 mbsf becomes supersaturated and gas hydrate forms. This depth is the upper edge of the zone with dispersed gas hydrate in sediments is described in the ODP report for site 997 [51]. As more methane is produced with time, more gas hydrate is formed. As demonstrated in the time-lapse video (Video S1), the gas hydrate produced by in situ organic matter degradation is less than 1% of the pore volume. With the continuation of deposition, gas hydrate is buried down as sedimentation proceeds and starts to dissociate when it crosses the base of the GHSZ depth, and after saturating pore According to the ODP report for the studied site [51], Blake Ridge contains gas hydrate, and below the GHSZ, there is an accumulation of free gas [16,57], which is produced via microbial methanogenesis [22]. According to the salinity profile, we considered up to 8% of pore volume could be occupied by gas hydrate and free gas. As shown in Figure 4, based on the model, methane produced by methanogenesis are in three phases, dissolve in pore water, as gas hydrate, and free gas. In the beginning, there is no methane in the system. By consuming all sulfate in depth, methane starts to be produced by methanogenesis. Over time, more methane is generated, and the concentration of methane in pore water increases, so it diffuses to the upper part of the column. After around three million years, when about 350 m of sediment accumulated, methane in pores below about 180 mbsf becomes supersaturated and gas hydrate forms. This depth is the upper edge of the zone with dispersed gas hydrate in sediments is described in the ODP report for site 997 [51]. As more methane is produced with time, more gas hydrate is formed. As demonstrated in the time-lapse video (Video S1), the gas hydrate produced by in situ organic matter degradation is less than 1% of the pore volume. With the continuation of deposition, gas hydrate is buried down as sedimentation proceeds and starts to dissociate when it crosses the base of the GHSZ depth, and after saturating pore water, free gas (methane bubbles) forms. On the other hand, due to the upward migration tendency of free gas by advection and diffusion, gas produced in a deeper part of the model starts to move up and enter the GHSZ to form more gas hydrate. Over time, more gas hydrate is formed and dissociated, forming a gas hydrate-rich level above and a gas-rich level below the base of the GHSZ. This contrasting horizon can be recognized in seismic sections as a bottom simulating reflections [40], which is a distinct feature in the Blake Ridge at the site 997 [47]. With time, the gas hydrate level becomes thicker and inhibits further upward gas migration, so gas accumulates below the GHSZ.

Methane Phase Behavior
Additionally, at the end of the modeled period (after 98.5 Myr), there is a slight decrease in pore volume occupied by gas hydrate because of a decrease in POC input (Figure 4 and time-lapse video in Video S1).

Carbon Sink
The degradation of POC, the only carbon input in the model, is transformed into methane (dissolved in pore water and as gas hydrate and free gas) and DIC, which is kept dissolved in pore water or reacts with calcium and magnesium to form carbonates [58,59]. Over time, DIC can diffuse to seawater from the top (seafloor); it can be removed from the bottom (600 mbsf) by advection (burial) or diffusion, the latter depending on the concentration profile. Dissolved methane can be removed from the bottom of the system by advection and removed or added to the system from the bottom by diffusion. Free gas can enter the system by advection and be removed or entered from the bottom by diffusion. Magnesium carbonate and calcium carbonate are in the solid phase and can be removed from the bottom by advection. There is no methane hydrate in the bottom of the modeled column because it is outside of the GHSZ (<450 mbsf). However, hydrate can be dissociated below the GHSZ, releasing methane. Time-lapse video for carbon balance is presented as Video S1. Figure 5 shows the quantity and distribution of consumed carbon per unit of the modeled column length after the last time step (i.e., 100 Myr). According to the model, about 17% of POC input is labile and consumed/mineralized in the system. From this 17%, carbon mass balance shows that DIC is the main carbon sink within the sediment (55.27%), and generated methane (dissolved, gas hydrate, and free gas) is the second carbon sink (23.13%). Methane hydrate and free gas host a minor part (0.6% and 0.9%, respectively) of the system's total carbon. However, as shown in Figure 3, the quantity of free gas trapped below the base of the GHSZ is significantly more than the value of 0.9%, the amount of free gas produced by in situ POC degradation, which does not include the free gas entered to the modeled column from below. water, free gas (methane bubbles) forms. On the other hand, due to the upward migration tendency of free gas by advection and diffusion, gas produced in a deeper part of the model starts to move up and enter the GHSZ to form more gas hydrate. Over time, more gas hydrate is formed and dissociated, forming a gas hydrate-rich level above and a gasrich level below the base of the GHSZ. This contrasting horizon can be recognized in seismic sections as a bottom simulating reflections [40], which is a distinct feature in the Blake Ridge at the site 997 [47]. With time, the gas hydrate level becomes thicker and inhibits further upward gas migration, so gas accumulates below the GHSZ. Additionally, at the end of the modeled period (after 98.5 Myr), there is a slight decrease in pore volume occupied by gas hydrate because of a decrease in POC input (Figure 4 and time-lapse video in Video S1).

Carbon Sink
The degradation of POC, the only carbon input in the model, is transformed into methane (dissolved in pore water and as gas hydrate and free gas) and DIC, which is kept dissolved in pore water or reacts with calcium and magnesium to form carbonates [58,59]. Over time, DIC can diffuse to seawater from the top (seafloor); it can be removed from the bottom (600 mbsf) by advection (burial) or diffusion, the latter depending on the concentration profile. Dissolved methane can be removed from the bottom of the system by advection and removed or added to the system from the bottom by diffusion. Free gas can enter the system by advection and be removed or entered from the bottom by diffusion. Magnesium carbonate and calcium carbonate are in the solid phase and can be removed from the bottom by advection. There is no methane hydrate in the bottom of the modeled column because it is outside of the GHSZ (< 450 mbsf). However, hydrate can be dissociated below the GHSZ, releasing methane. Time-lapse video for carbon balance is presented as Video S1. Figure 5 shows the quantity and distribution of consumed carbon per unit of the modeled column length after the last time step (i.e., 100 Myr). According to the model, about 17% of POC input is labile and consumed/mineralized in the system. From this 17%, carbon mass balance shows that DIC is the main carbon sink within the sediment (55.27%), and generated methane (dissolved, gas hydrate, and free gas) is the second carbon sink (23.13%). Methane hydrate and free gas host a minor part (0.6% and 0.9%, respectively) of the system's total carbon. However, as shown in Figure 3, the quantity of free gas trapped below the base of the GHSZ is significantly more than the value of 0.9%, the amount of free gas produced by in situ POC degradation, which does not include the free gas entered to the modeled column from below.   by in situ organic matter degradation, respectively. Subscribes "top" and "bottom" stand the amount of substance that leaves the modeled column from top and bottom, respectively. Subscribe "labile" indicates the amount of converted (labile, mineralized) organic matter.

Sedimentation Rate and Initial Age of POC
To find the effect of sedimentation rate in biogeochemical processes, site 997 with constant POC at the top of the column was modeled. Since a higher sedimentation rate is equivalent to younger organic matter age [6,9,15], we considered three different initial ages with a linear function of sedimentation rate. The three sets are: (a) w = 0.01 cm/yr and a = 8 × 10 5 yr (oldest), (b) w = 0.02 cm/yr and a = 4 × 10 5 yr (middle-aged), and (c) w = 0.04 cm/yr and a = 2 × 10 5 yr (youngest)). The results of the last time step of each set are plotted in Figure 6. By decreasing POC initial age (increasing the sedimentation rate), the organic matter becomes more reactive, and the POC degradation rate increases (according to the results, the rate of youngest POC is almost 3.73 and 1.93 times higher than the rate of the oldest and middle-aged one, respectively) results in increasing sulfate reduction and methanogenesis rates (Figure 7), and consequently, shallowing the SMT and increasing methane production ( Figure 6). As methane generation increases, pore water is saturated at shallower depths, hydrate formation commences in shallower depths, and the thickness of hydrate (the supersaturation zone of methane) increases (276, 330, and 366 m for the oldest, middle-aged, and youngest, respectively), which is in general agreement with other studies [9,60]. Nevertheless, the gas hydrate average concentration decreases with an increase in sedimentation rate. Upward gas migration velocity in the lowest sedimentation rate is the highest; thus, more gas reaches GHSZ, and consequently, more hydrate is formed at the same time. On the other hand, the gas hydrate in the system with the lowest velocity is buried slower, and it leaves the GHSZ more slowly; therefore, it dissociates with a slower rate compared to the other two modeled systems. As a result, the amount of gas accumulated below the GHSZ is not proportional to the hydrate amount.

Coexistence of Free Gas and Gas Hydrate in GHSZ
Here we propose an idea to explain the coexistence of free gas and gas hydrate in GHSZ from the kinetic perspective. The presence of materials in pore water and sediments with the ability to retard gas hydrate formation and suppress the gas hydrate growth can decline the gas hydrate formation rate (such as polysaccharides [61,62] and polyvinylcaprolactam [63] as natural and artificial kinetic hydrate inhibitors, respectively). To investigate the effect of gas hydrate formation kinetic in pores, after about 37 million years, we assumed the kinetic constant of gas hydrate decreases to 0.1, 0.01, and 0.001 of the initial value. The model was run for over 50 million years. By reducing kinetic constant, the rate of hydrate formation decreases, and since the sedimentation rate is constant, hydrate dissociates faster compared to the formation, so the pore volume occupied by gas hydrate diminishes, and the excess (supersaturation) of methane is released as bubbles (see time 38 Myr in time-lapse video in Video S2). On the other hand, as shown in Figure 8, by decreasing the gas hydrate kinetic constant, free gas below the GHSZ can migrate upwards to shallower depths. According to the model, free gas in GHSZ causes pore water to be saturated at shallower depths, increasing gas hydrate thickness and shallowing the SMT. Although the rate of AOM increases (Figure 9), it is not enough to consume all the methane, and some part of free gas can vent out the seafloor. The efficiency of AOM barrier is 97.75, 38.78, and 28.22% for the gas hydrate kinetic constants of 0.1, 0.01, and 0.001 K GH , respectively. This indicates that, by increasing the amount of free gas in top layers, AOM cannot consume all methane and leads to methane efflux to the sea-water interface. There is indication that dissolved organic carbon [64] and potentially some kinetic inhibitors may be concentrated in the sediments over time. As a suggestion, we recommend to investigate the existence of the kinetic gas hydrate inhibitors in future expeditions, especially at areas with seafloor methane venting.

Conclusions
We used a dynamic model to simulate the biogeochemical process in marine sediments related to the mineralization of organic matter, including methane formation and consumption. Compared with other studies, our model predicts dissolved species' temporal and spatial concentration in pore water, considering the seawater concentration as the top boundary condition without any subsurface information. Prediction of methane generation and consumption, SO 4 2− , Ca 2+ , Mg 2+ , and DIC with reliable accuracy is an important achievement of our model. The model predicts that in site 997, the most gas hydrate accumulated in the sediment is formed from free gas generated by organic matter mineralization below the GHSZ. Based on our model, a massive volume of free gas below the GHSZ could migrate upwards to shallow depth if the gas hydrate formation kinetic is adequately slow.
Although our model fitted the data well by using the seafloor data as the only input, the model can be further improved by considering the energy balances to predict the effect of global warming on the methane emission in the oceans and hydrate stability.
As a future study, the model could also be used as a new tool to estimate the total carbon budget and distribution in the sediments in areas with limited data. The prediction of the methane accumulation in sediments is another clear implication of the model for exploring energy resources and understanding the past and future of Earth's climate. Moreover, the authors are going to investigate the effect of temperature, salinity, and eutrophication on biogeochemical processes. In the following, a detailed description of the model is presented in Tables A1-A4. In this study, the porosity is considered to decay exponentially with depth due to steady-state compaction [50], the porosity at depth zero (∅ 0 ), the porosity at infinite depth (∅ f ), and the attenuation coefficient (p) which are calculated based on the available data in the site. Decrease of the porosity because of the hydrate formation is considered in the gas hydrate stability zone by distracting the volume of hydrate from the pore volume. Burial velocity of pore water and solids are calculated based on steady-state compaction [50] and the available data is used for the sedimentation rate (w f ). The upward velocity can be in the range of 0.005 to 30 cm yr −1 [4]. In this study, we assumed the upward velocity of gas as 0.02 cm yr −1 . Molecular diffusion coefficients (D M ) are calculated in different temperature and salinity according to depth while logarithmic equation is used to consider the effect of the porosity [52].
The rate laws and expressions used in the modeling are described in Tables A2 and A3. Modified Middelburg's age-dependent model is used as POC degradation rate, where K C is inhibition constant to decrease the rate of reaction with increase in metabolites (CH 4 and DIC), and a 0 is the initial age of the organic matter [8,9]. The rate of sulfate reduction and methane formation are related to organic matter degradation rate, where K SO 4 is a Monod constant to decrease the sulfate reduction rate at low concentration of sulfate (K SO 4 = 1 mM) [9]. A simple second-order kinetics is used for the anaerobic oxidation of methane (AOM) [39], where K AOM is a kinetic constant. The calcium carbonate precipitation and dissolution rates are considered to depend linearly on the saturation state, where K S,CaCO 3 is the thermodynamic equilibrium constant modified by the pressure [39,65,66]. The magnesium carbonate precipitation and dissolution rates are assumed to depend linearly on the saturation, same as the calcium carbonate, where K S,MgCO 3 stands for thermodynamic equilibrium constant modified by the pressure [65,67]. In this model, whenever the concentration of dissolved methane surpassed the methane solubility (C CH 4 ,sat ), hydrate is precipitated if it is in the hydrate stability zone; otherwise, free gas is released. The methane solubility is calculated by modified Pitzer-approach [68] and Duan-approach [69,70], in and out of hydrate stability zone, respectively. Hydrates and free gasses were allowed to dissolve in pore water if the concentration of methane is less than saturation. The kinetic parameters of the model were constrained by fitting the model to the porewater profiles (Table A4).

Parameter Equation
Burial velocity of solids w = w f (1−∅f )  Table A2. The rate laws used in the modeling.

Rate Expression
POC degradation R POC =  Table A3. Rate expression used in the differential equations.

Species Rates
Particulate   Figure A1. Spatial and temporal results of the model for the simulated ODP site 997 [49]. First row: Concentration profile of (a) Dissolved methane (CH4); (b) Methane in the gas phase (CH4 (g)); (c) Methane in hydrate phase (GH). Second row: Simulated rate of (d) Methanogenesis (MP); (e) Free gas formation (CH4 (G)), and (f) Gas hydrate formation (GH). Last row presented in shorter depth (to magnify the rate changes): (g) Concentration profile of sulfate (SO4 2− ); simulated rate of (h) Sulfate reduction (SR), and (i) Anaerobic oxidation of methane (AOM). The black area refers to the place where no sedimentation occurred.