A two-dimensional model of the methane cycle in a sedimentary accretionary wedge

A two-dimensional model of sediment column geophysics and geochemistry has been adapted to the problem of an accretionary wedge formation, patterned after the margin of the Juan de Fuca plate as it subducts under the North American plate. Much of the model description is given in a companion paper about the application of the model to an idealized passive margin setting; here we build on that formulation to simulate the impact of the sediment deformation, as it approaches the subduction zone, on the methane cycle. The active margin configuration of the model shares sensitivities with the passive margin configuration, in that sensitivities to organic carbon deposition and respiration kinetics, and to vertical bubble transport and redissolution in the sediment, are stronger than the sensitivity to ocean temperature. The active margin simulation shows a complex sensitivity of hydrate inventory to plate subduction velocity, with results depending strongly on the geothermal heat flux. In low heat-flux conditions, the model produces a larger inventory of hydrate per meter of coastline in the passive margin than active margin configurations. However, the local hydrate concentrations, as pore volume saturation, are higher in the active setting than in the passive, as generally observed in the field.


Introduction
An accretionary wedge sediment complex is an example of a "structural" hydrate deposit, in which tectonics and gas flow play obvious roles in controlling the abundance and distribution of methane hydrate, as opposed to the "stratigraphic" deposits (Milkov, 2004), in which sediment accumulates into depositional layers.In an accretionary wedge complex, the sediment is actively deformed by compression associated with the scrape-off of the sediment complex from the underlying subducting oceanic crust (Carson et al., 1990).
Structural deposits comprise a smaller fraction of the hydrate inventory globally (Boswell and Collett, 2011), but form hydrate in higher concentrations, even forming massive hydrate blocks.Hydrate is found closer to the sea floor than in the stratigraphic deposits (Torres et al., 2002), which tend to concentrate hydrate at the base of the stability zone, often hundreds of meters below the sea floor.Structural deposits might therefore release carbon into the ocean or even methane into the atmosphere more quickly than the time constant would be for a response from the stratigraphic deposits.
Models of methane hydrate cycling have been applied to active margin settings by driving them with increased upward pore fluid flow (Buffett and Archer, 2004;Chatterjee et al., 2011;Luff and Wallmann, 2003).Here we attempt to internalize the pore fluid flow into a methane hydrate sedimentary model by expanding the model domain to bedrock and into a second dimension, perpendicular to the motion of the underlying plate sliding into the subduction zone.In particular, the model is applied to the case of the Cascadia margin (Spence et al., 2000).

Overview of the passive margin configuration
Called SpongeBOB, the numerical model is described in a companion paper (Archer et al., 2012) as it was formulated for a passive continental margin setting.The model formulation as described in that paper will be summarized here before we show details of the additional model formulation Published by Copernicus Publications on behalf of the European Geosciences Union.

D. E. Archer and B. A. Buffett: A two-dimensional model of the methane cycle
required for the accretionary wedge setting.The model is formulated on a two-dimensional grid, onshore/offshore in the lateral dimension and with a stretching "sigma" grid in the vertical.The model is intended to span the continental margin from the continent to the abyss, over geologic time scales of 10 7 -10 8 yr.

Sediment transport
A sediment transport scheme distributes material to the sea floor.Continental material originates at the left-hand (continental) side of the domain, and is transported and sedimented according to the sinking velocities of the various grain sizes, and the water depth (which prevents sedimentation if it is too shallow).Another fraction of the sedimenting material is called "pelagic", and it deposits uniformly throughout the domain, if the water is deep enough.When the slope of the sea floor exceeds a critical value set at 4-6 % grade, sediment resuspends and is distributed downslope.This material is assumed to carried to the abyssal floor in turbidity currents that only allow resedimentation to begin when the slope of the sea floor decreases offshore to less than 1 % grade (Meiburg and Kneller, 2010).The sedimentation scheme produces a continental shelf, a well-defined shelf break, and a continental slope.The various parameters of the sedimentation were tuned in order to reproduce the envelope of sediment: the shapes of the sea floor and the depth to bedrock.

Isostasy
Bedrock in the model floats isostatically, balancing the load from the crust and sediment against that of a hypothetical displaced mantle fluid.The buoyancy of the crust is affected by cooling of the upper mantle in a thermal boundary layer that thickens with the square root of time, allowing subsidence with increasing crustal age.The elevation of the crust relaxes toward the equilibrium value on an isostatic rebound timescale of 10 4 yr.The passive margin model has no representation of crustal rigidity except for a numerical smoothing operation, which has a spatial range of 10-20 km.

Organic carbon and methane
The particulate organic carbon (POC) content of the continentally derived sedimenting material is specified in the model when it hits the sea floor, as a function of water depth.In the passive margin simulation, sea level changes were imposed on the simulation, along with a correlated time-varying oxygen "state" of the ocean, which drives changes in POC concentration and the chemistry of the organic matter, most notably its H / C ratio.The active margin simulation reaches steady state in only 10 Myr, which is shorter than the 200 Myr duration of the passive margin simulations, so geological sea level and ocean oxygenation changes are not imposed on the simulations here.Instead, the relative sea level for these simulations was varied by ±20 m on a cycle time of 1 Myr, rep-resenting local tectonic uplift and subsidence driven by the forces confronting the crust in this "turbulent" part of crustal geophysics.The impact of this stipulation can be assessed in a simulation called No Sealevel with time-invariant sea level.
Biologically and thermally driven chemical reactions produce dissolved methane, CH 4 , which interacts with the gas and hydrate phases depending on temperature and pressure conditions.Respiration of POC first consumes pore water sulfate, SO 2− 4 , until it is depleted, which occurs relatively shallow in the sediment column.Bacterial respiration of POC then produces CH 4 and CO 2 .The maximum efficiency of CH 4 production from the organic carbon is set by redox balance according to the H / C ratio of the POC.Porewater δ 13 C data constrain the methanogenesis to be about 50 % of respiration, which is a bit lower than the maximum set by the redox constraint, as if some of the molecular hydrogen intermediary reacts with oxidized mineral phases rather than dissolved CO 2 to produce methane.
As temperatures warm further, exceeding about 60 • C, petroleum is produced, if the H / C ratio in the POC exceeds a value of 1 (Hunt, 1995).Petrogenesis draws the H / C ratio down toward a value of 1.A fraction of this petroleum (10 %) is assumed to migrate upward with a velocity of 1 m per thousand years.If it reaches the biological zone (temperature less than about 50 • C), it can be respired, producing CH 4 and CO 2 , similarly to respiration of POC.
Thermal methanogenesis begins at about 150 • C, producing CH 4 , dissolved CO 2 species collectively called dissolved inorganic carbon (DIC), and dissolved organic carbon (DOC).Ultimately, the CH 4 production is limited by hydrogen in the POC, with H / C approaching 0 in the hottest sediments.Thermogenic DOC production is produced in a stoichiometry of CH 2 O (e.g.acetate).DOC is also released by "sloppy feeding" in the respiration zone, and it is consumed in the respiration zone (producing DIC and methane if SO 2− 4 is not available) with the same rate constant as applied to migrated petroleum.SO 2− 4 and CH 4 are also consumed by anaerobic oxidation of methane (AOM).
The CH 4 concentration is compared with the solubility of bubbles and hydrate, and allowed to form those phases if it exceeds supersaturation.Hydrate is stationary within the sediment but bubbles migrate, following an ad-hoc parameterization that redistributes CH 4 in bubbles into overlying grid cells.This parameterization was needed to prevent buildup of excessive bubble volumes throughout the deep sediment column, and, although the details of how the transport occurs are sketchy and at any rate unresolved in the model, it seems clear that, in the absence of evaporates or permafrost, most methane gas produced in the sediment column does manage to escape the sediment column eventually (Hunt, 1995).
When rising gas encounters undersaturated conditions in the model, it redissolves following a simple exponential function of height in the sediment column.The methane inventory was found to be extremely sensitive to the scale height in the parameterization, and we will show similar sensitivity studies for the active-margin model with similar results.

Reference case and variants
To understand and document how the model works, we show results from a suite of model sensitivity runs summarized in Table 1.Details of these scenarios will be explained as the new components of the model relevant to the active margin setting are described.In contrast to the passive margin simulations, the grid resolution is the same for all of the active margin simulations presented here.

Deformation of the sediment column
In the active model configuration, in addition to the processes in the passive margin model, the model grid in the horizontal dimension is manipulated to simulate the uniform compaction and thickening of the sediment column by lateral compression.The x coordinate values of the grid cells are carried laterally by the moving crust, and the spacing between the grid points decreases, resulting in uniform vertical thickening of the sediment column.
The velocity of the incoming sediment column in the offshore edge of the model domain (the right) is the crustal velocity, specified as a parameter of the model scenario (for which there are sensitivity runs Plate 100, 80, and 20).The base plate velocity is taken to be 40 mm yr −1 , from the subduction rate of the Juan de Fuca plate.The Juan de Fuca Ridge is neither orthogonal to the direction of plate motion, nor is it geographically stationary, but rather is moving slowly toward the trench.The model formulation as presented here simplifies this geometry into two dimensions by equating the spreading and subduction velocities, maintaining a constant distance between the ridge and the subduction zone throughout the model simulation.This simplification affects the plate velocity and also the amount of sediment that enters the wedge through time, and its impact can be assessed from the sensitivity to plate velocity and to sedimentation (simulation Pelagic).
At the onshore end of the domain, the sediment column velocity is held fixed at a value 10 times lower than the incoming sediment column velocity, representing a sediment column nearly stopped by collision into the other plate.The rate is specified at an extrapolated x location at which the sea floor would rise above sea level.The velocities (u) of the grid points offshore of this are determined by where K deform is a dimensionless deformation constant that controls the horizontal extent of the wedge relative to the thickness of the incoming sediments.The velocities are negative because they flow in the direction of decreasing x (Fig. 1a).The horizontal extent of the deformation zone, in this simulation about 150 km, is determined by the value chosen for K deform (for which there is a sensitivity run called Wide Def ).The x coordinate (location in physical space) of each horizontal grid point is updated each time step according to its calculated velocity.As the sediment column slides landward, the state of the underlying ocean crust at the grid point is recalculated based on its new location, including in particular the flexure effect of the subducting plate in the trench.The sediment column rides but also isostatically steers the pathway of the subsiding ocean crust as it thickens.
As grid points move across the x = 0 origin, or as the sediment column begins to outcrop from the ocean, they are dropped from the model domain, and a new grid point is created on the far right-hand side of the domain.These sediment columns are initialized with a computationally required minimum 1 m thickness per each of 15 grid cells, a negligible fraction of the ultimate sediment wedge.The 400 km mark from the extrapolated coastline is taken to be the spreading center where ocean crust is created and begins to accumulate sediment.The grid points propagate through the domain, emerging on the right as new crust is formed at the spreading center, and accumulating sediment and deforming as they converge toward the left.The eventual model steady state has no stationary points in it, but it manages to reach a moving stationary solution.
By conservation of volume, the thickness of the column z sedcol increases as x decreases.The model is formulated in vertical columns, which requires that the sediment column deforms strictly vertically.The increasing vertical has the effect of increasing the excess pressure in the fluid phase, driving an expulsion fluid flow (Yuan et al., 1994).Real sediment column compression is generally focused on diagonal faults in the column, with a block of sediment from one side of the fault over-riding the other by sliding upward along the fault.The faults can dip onshore (normal) or offshore (abnormal), perhaps depending on the frictional state of the contact with bedrock (Davis et al., 1983).
Fortunately, previous models of sediment accretion have found a strictly vertical formulation to be an acceptable approximation, for modeling heat flow (Wang et al., 1993) and mineral closure ages (Batt et al., 2001).The effect of slip motion along diagonal faults dipping in either direction would be to displace material laterally.A parcel of over-riding sediment over a fault dipping offshore, for example, would be moving shoreward somewhat as it rides the wedge upward toward shallow waters through the domain of the model.A parcel at the top of the undeformed incoming sediment column might progress into the wedge a bit before a contemporaneous parcel from the bottom of the incoming column.But the lateral displacements ought to be limited by typical fault geometry to be not much larger than the thickness of www.biogeosciences.net/9/3323/2012/Biogeosciences, 9, 3323-3336, 2012 Broad Slope Critical seafloor slope before sliding of 2 % instead of default 6 %.

Bumpy
Heterogeneous sediment deformation constant imposed as 100 km variations of 40 % in the deformation constant.
Wide Def Decreased sediment deformation constant, spreading the deformation zone from 8 × 10 −2 in the Base scenario to a value of 4 × 10 −2 .The change diminishes the need for erosion to maintain a critical seafloor slope.
No Erode Erosion disabled.
No Thermogen Thermogenic methane production disabled.

No Chan
No vertical low-permeability chimneys.
No Sealevel Sea level oscillations of ±20 m on 1 Myr time cycle in Base simulation disabled.the sediment column, 5-10 km or so.This is a small displacement relative to the overall width of the wedge, which is about 100 km.Therefore, we expect the effects of this mode of motion to be relatively small.Ultimately, the inventory of sediment in the model domain overall is determined by the balance of the sources and sinks, between sediment deposition from the adjacent continent and lateral sediment advection out of the model domain.Choosing the factor by which to impede the outgoing sediment column velocity (the factor of 10) essentially sets the model domain; a value of 10 achieves a solution in which the sediment column is close to outcropping at the sea surface, encompassing the entire hydrate stability zone but missing the complexities of sediment transport and erosion that produce the continental shelf, and erosion on land (Fig. 2).
The model as described so far produces a smooth continental slope, but the sea floor in real accretionary zones is ridged, as blocks ride over each other.We attempt to simulate the effect of topography in the wedge by varying the deformation constant laterally, by ±40 % on a wavelength of 100 km, in a simulation called Bumpy.The zones of high and low deformability travel with the material through the domain, and new grid points are initialized by extrapolation, continuing the original wave into the new incoming model domain (Fig. 2).

Crustal bending
Near the subduction zone, an oceanic plate is pressed downward by the load of the subducted lithosphere on the other side.The lithosphere deforms elastically, with the flexural rigidity determining the bending in response to a given torque.The situation is analogous to a floating dock with a person ready to dive in, standing at the edge.There is a zone next to the diver where the lateral cohesion of the crust pulls it downward, and then inshore of this, the dock rises out of the water somewhat, due to the requirement for overall isostatic equilibrium, and to the stiffness of the dock.This situation is treated in the model based on an analytical solution to the case of a single point load at the subduction zone (Turcotte and Schubert, 1982).The differential equation is

Biogeosciences
and its solution with application of appropriate boundary conditions where There are three tunable parameters in this formulation (Fig. 1b).One, z b , corresponds to the height above isostasy of the forearc bulge, for which we use 100 m.Two are horizontal space scales.One, x 0 , is from the edge of the plate (where the guy is standing, x = 0) to the boundary between the trench and the bulge (defined as the local isostatic equilibrium z = 0 line).The other, x b , is the coordinate at the peak of the bulge.We use 125 and 175 km, respectively.These horizontal scales have a huge impact on the eventual depth of the trench.The isostatic load of the sediment column, described next, greatly amplifies the eventual depth of bedrock in the trench.A "torque pulldown" of about 2000 m at the left-hand side of the domain ultimately results in 10 km of eventual depth to bedrock (Fig. 2).

Isostasy
The displacement of the crust near the subduction zone is implemented in the model as a deviation from the crust unloaded by sediment as described in Sect.3.3.The correct way to do this calculation would be to solve the differential Eq. ( 1) using the distributed load of the sediment column, allowing the load full interplay with the rigidity of the crust.SpongeBOB is formulated as a simpler approximation of this.The mass load affects the isostatic equilibrium value locally, without regard for the springiness of the crust, in order to benefit from the convenient analytical solution to the x = 0 point load case.In this way, the mass of the sediment column still has an impact on the elevation of the crust, so that sedimentation can drive subsidence.

Sea level and time-dependent forcing
Because of the shorter simulation time (10 Myr), the 150 Myr sea level cycle, a fundamental driver to the passive margin simulations (Archer et al., 2012), was neglected.We do however incorporate a ±20 m sea level oscillation on a time scale of 1 Myr to simulate the turbulence that a tectonic regime must be experiencing.The sensitivity to this forcing is gauged by a simulation called No Sealevel.

Sediment erosion and landslides
The solution of the horizontal compression and deformation scheme, by itself, would tend toward ever-increasing sea floor slope in the shoreward direction (Fig. 3).This tendency is balanced in the model by slope erosion and landslides.The grade of the sea floor is limited to a critical value (6 % in the Base simulation) by two mechanisms.One is during sediment deposition; if the sea floor slope is supercritical, the material that would have sedimented is instead added to a resuspended pool and advected offshore.
The other mechanism is an erosional term that is triggered when the sea floor slope exceeds critical.Sediment is removed from the top computational box by relaxation toward a value that would bring the grid point back to the critical slope with respect to its adjacent grid point, with a relaxation time constant of 10 −5 yr −1 , sufficiently fast to hold the sea floor close to the critical value even while the sediment column is steepening by deformation, but slow enough to be kind to the numerics of the model.Solid material and pore fluid are advected upward through the computational grid in the interior of the sediment column, allowing the expandable "sigma" grid to contract as the sediment column erodes.The resuspended material is assumed to be incorporated into a turbidity current that travels down slope without deposition of sediment until the sea floor slope is less than 1 %.Beginning at this point, the resuspended sediment de-posits on the sea floor following the same criteria for deposition as used by the primary depositing continental material, with larger size classes falling out faster than small particles.Snapshots of sediment accumulation/erosion rates and the fraction of redepositing material are shown in Fig. 4.
The eroding material conserves the chemical and grainsize characteristics, but the redeposition of POC is fractionated by the size separation of the redepositing material.The POC is distributed, in the model as well as in observations (Mayer, 1994), according to the surface area of the sediment, such that smaller size classes have a higher POC content by weight, due to their larger surface-to-volume(mass) ratio.The size fractionation of the POC can move the zone of highest POC content offshore; in particular this is evident in the Broad simulation, which has a shallow critical slope angle and hence a lot of sediment redeposition (Fig. 5).

Pore fluid flow
Pore fluid flow in the model is driven by the accumulating mass of the sediment column, and governed by Darcy's law and a sediment permeability that depends on the grain size and the porosity of the local sediment.The derivation of excess pressure from the porosity and the numerical advection scheme, including flow limiter, are described in the companion paper (Archer et al., 2012).Most of the fluid flow from the real ocean sediment column appears to make its way through high-permeability channels and pathways rather than flowing homogeneously through the bulk sediment column.The SpongeBOB model is too coarsely gridded to resolve faults and sandy turbidites in detail, but the overall impact of flow heterogeneity on the geochemistry of the sediment column is mimicked by creating vertical channels of high permeability within the SpongeBOB grid.The channels are placed every 5 model grid points, and their permeability is enhanced by a factor of 10 relative to the background grid points.The grid points are found to focus the upward flow, and result in significant changes in the pore water chemistry both within the channels and also in the background cells.For these active margin simulations, the permeable channels are carried horizontally with the computational grid, following the sediment material.Vertical flow velocities, relative to the sediment grains (Darcy flow, w Darcy ) and relative to the sea floor (w seafloor ), are shown in Fig. 6.

Heat flow
The heat flow results from the sediment surface are shown in Fig. 7.The model captures the general trend of lower heat fluxes onshore, and by about the same magnitude of drawdown.The model fluxes are slightly lower than the measurements from Hyndman and Wang (1993), driven by a model geothermal heat flow of 120 mW m −2 , compared to the offshore heat flux of 120-130 mW m −2 in the data.The simulations all reproduce the trend of decreasing heat flow with sediment column thickening observed in the data.The diffusive heat flow values are impacted by the permeable vertical channels, which appear as spikes of high heat flux.
The one curious model result comes from the Broad scenario, in which a shallow critical slope angle drives high rates of sediment erosion and redeposition (Fig. 7b, dashed lines).
In this scenario, removal of sediment by erosion increases the heat flux landward, counteracting the decrease landward seen in all the other simulations (and the data).

Organic carbon and methane
The model respiration and thermal degradation kinetics are the same as in the companion paper (Archer et al., 2012).
The production rates of methane from respiration and thermal degradation are shown in Fig. 8. Methane concentrations from the simulations are shown in Fig. 9. Bubbles are shown in Fig. 10, and hydrate in Fig. 11.The highest hydrate saturations tend to be found at the toe of the wedge, similarly to observations from Cascadia (Malinverno et al., 2008).The simulations are similar to the eye, but variations can be seen in the envelope of the sediment column, for example, due to plate velocity (Plate 100, 80,and 20), sediment column deformation (Wide Def and Bumpy), and sediment transport and reposition (Pelagic and Broad Slope).Methanogenesis rates are affected by the labile fraction of POC (Bio 10 % and Bio 100 %), impacting the distribution of dissolved methane, bubbles, and hydrate.
Stable carbon isotopes provide a diagnostic window into the dynamics of the sedimentary carbon cycle.Concentrations of dissolved inorganic carbon (DIC) are shown in Fig. 12, and profiles of the δ 13 C of dissolved methane and dissolved inorganic carbon (DIC) from the Base case are compared with the measurements of Pohlman et al. (2009) in Fig. 13, and results from the other scenarios are plotted in Figs. 14 and  model geophysical scenarios do not have a major impact on the isotopic compositions.
In broad brush the model captures the general isotopic compositions as measured, but the field data show a systematic dependence of the isotopic compositions that is not found in the model results.The cores were drilled in a transect from the offshore toe of the wedge, in water depth of about 2200 m, to water depths of about 1000 m in the inshore direction, about the middle of the wedge.The isotopic compositions of DIC and CH 4 are both systematically heavier in the inshore core than they are near the toe of the wedge.The differences apply in parallel to both tracers, as though the isotopic composition of methane in surface sediments were controlled by constant fractionation offset from the DIC parent material.However, the methane isotopic composition in near-surface sediment can also be impacted by fractionation during anaerobic oxidation of methane by sulfate (AOM by SO 2− 4 ), and methane isotopic compositions in the model tend to deteriorate as the methane becomes depleted near the sediment surface.The isotopic composition of DIC at depth in the model sediment is governed not by the rate of methane formation, but rather by the proportion of the reacting organic carbon that produces CH 4 versus DIC.Methanogenesis from organic matter ultimately acts as a source for DIC, of an isotopic composition that depends on the relative proportions of DIC and CH 4 produced.For example, if the source POC has a carbon isotopic composition of −25 ‰, and the fractionation on CO 2 reduction is −60 ‰, then a conversion efficiency of POC to CH 4 of 25 % would explain the isotopic compositions of both tracers at the toe of the deformation zone, and an increase to 75 % methanogenesis efficiency would explain the measurements in the mid-wedge.The model does not suggest an explanation for why the methane production efficiency should vary in this way, however, as the isotopic compositions in the model are relatively invariant between the toe and the middle of the sediment wedge.

Model hydrate inventory sensitivities
Figure 16 shows the hydrate inventories of all the simulations, and the results are digested into specific model sensitivities, and compared with results from the passive margin simulation, in Fig. 17.
The sensitivity of the model to ocean temperature is similar between the passive and active models (Fig. 17a).Both configurations of SpongeBOB are less temperature sensitive than the Davie and Buffett (2001) model as deployed globally by Buffett and Archer (2004).
Many of the other model sensitivities can be rationalized as resulting from the differing sediment dynamics within the hydrate stability zone in the two model configurations.The residence time of sediment in the hydrate stability zone is much shorter in the active margin simulation, because sediment flows horizontally through the wedge on a time scale of only a few million years.In the passive margin setting, the lifetime of solid sediment in the stability zone is determined by sediment accumulation, resulting in a residence time an order of magnitude longer.
As a sediment parcel is carried laterally through the wedge, its depth below the sea floor increases as the sediment column overall stretches in the vertical, carrying it for example from shallow sediment depth in a deep water setting, within the hydrate stability zone, to relatively deeper depths below the sea floor in the thicker sediment column in the shallower water onshore setting.The sediment is thus carried through the same stages as a sediment parcel in the passive margin configuration, from within the stability zone to below the zone, except on a horizontal pathway rather than vertical sediment accumulation.As the parcel crosses the hydrate stability boundary, on either pathway, the methane solubility reaches a maximum, then decreases as the parcel moves deeper.The fluid tends to degas as the methane solubility decreases below the stability boundary.The horizontal motion thus tends to act as an accumulator for methane, additionally and operating more quickly than the analogous methane accumulator from sedimentation and burial.At the same time, the shorter residence time for solid sediment in the hydrate accumulation zone tends to limit the response of the steadystate hydrate inventory to many of the driving parameters.
The active margin configuration requires a higher minimum availability of organic carbon (labile POC fraction   in these results) in order to cross the threshold to forming methane hydrate (Fig. 17b).Presumably, the shorter lifetime of sediment in the active margin stability zone requires somewhat higher rate of methanogenesis in order to reach saturation in the shorter allotted time.This is the one example of heightened sensitivity in the active margin configuration of the model.The model sensitivity to the bubble redissolution scale height (Fig. 17c) is extremely muted in the active margin setting.Bubble redissolution as it migrates upward in the model is a recycling mechanism, strengthening the efficiency of the hydrate stability zone as a cold trap.A more efficient cold trap (shorter redissolution scale length) has a greater impact in the passive margin setting, because the sediment lifetime is longer so the trap has more time to operate and accumulate hydrate.
The active configuration is also less sensitive to geothermal heat flux (in Fig. 17d).Heat fluxes tend to be higher in active margin settings than passive, but the model results show that the active margin configuration is more able to tolerate higher heat fluxes than the passive configuration.In the passive margin, a high geothermal heat flux creates a thinner stability zone, which in the fullness of time on the passive margin allows the diffusive loss of methane sufficient to deplete the hydrate inventory.In the active margin, the vertical deformation tends to decrease the temperature gradient through the wedge, diminishing the impact of the basal heat flux on the near-surface temperature gradient and thickness of the stability zone.
The model sensitivity to the plate velocity (Fig. 17e) appears to be a convolution of several effects.Results are shown from two different values of the geothermal heat flux: one from the Base case of the passive margin setting; and the other from the active margin.Recall from Fig. 17d that the model is much more sensitive to the geothermal heat flux in the passive margin, while the active margin appears to buffer the hydrate inventory against depletion under conditions of high geothermal heat flow.This same effect explains most of the plate velocity sensitivity results in Fig. 17e.The order-ofmagnitude decrease in the residence time of sediment in the hydrate-bearing zone manifests itself as a sharp change in the hydrate inventory between the passive margin results (0 plate velocity) and the active margin results.However, the hydrate inventory is much less sensitive to the heat flow variations within the active margin model, which affect by a smaller but still seemingly significant order of factor of two.
Other insights into the functioning of the model can be gleaned by examination of a few more of the simulation results in Fig. 16.The vertical permeable channels lead to a minor increase in the overall abundance of methane hydrate (No Chan) in both passive and active margin configurations.Thermogenic methane production is found to be unimportant to methane hydrate in the active simulation (No Thermogen), a smaller sensitivity than the factor-of-two decrease in the passive margin results.Thermal methanogenesis alters the δ 13 C of methane in the deep sediment column, but it has limited impact on compositions in the shallower sediment column (Fig. 14).Most of the geophysical variants had only a small impact on the hydrate abundance, e.g. the heterogeneous deformation scenario Bump, the altered sedimentation scheme Pelagic, or the No Erode scenario, which prohibits sediment erosion.The simulation No Sealevel shows that the fluctuations in sea level that are imposed on the Base case had little effect on the hydrate inventory.

Conclusions
Active margin coastal settings tend to have high concentrations of methane hydrate in surface and near-surface sediments, leading to a conceptualization of intense methane delivery to surface sediments by pore fluid flow, driven by the deformation of the thick sediment wedge.Previous models of the dynamics of methane in sediments were applied to active margins by increasing the upward vertical pore fluid flow in the model (Buffett and Archer, 2004;Chatterjee et al., 2011;Luff and Wallmann, 2003).An increase in upward flow leads to an unambiguous increase in methane hydrate inventory in this type of model.
The SpongeBOB model internalizes the pore fluid flow, predicting the flow rates from sediment accumulation, deformation, and pore pressure dynamics, rather than imposing them as independent driving variables.As such it is difficult to isolate the effect of pore fluid upward flow, because the rate of flow is driven largely by the rate of sediment column convergence, which also changes other factors such as the residence time of solid material in the hydrate-bearing zone.The model tends to produce higher pore-volume saturations (concentrations) of methane hydrate in the active margin configuration than it does in the passive, but the active margin does not necessarily contain more hydrate per meter of coastline (total inventory) than the passive margin configuration.The horizontal transit of sediment through the hydrate-bearing zone appears to buffer the hydrate inventory, making it generally less sensitive to various driving parameters such as the geothermal heat flux and the efficiency of the stability zone as a methane trap (the redissolution scale height parameter).
Heat 90, 150, 180 mW Geothermal heat flow in units of mW m −2 , relative to a base case of 120 mW m −2 .Pl xx Ht 70 Effect of plate subduction velocity (xx mm/yr) under conditions of low geothermal heat flux (70 mW m −2 ) for comparison with the passive margin configuration of the model.
Fig. 1.(a) The sediment column velocity (negative meaning from right to left) for the Base scenario.(b) The solution to the plate flexure/isostasy balance near the subducting margin.

Fig. 2 .
Fig. 2. Grid, temperature, and isostasy results for (a) Base, (b) Bumpy, (c) Broad Slope, and (d) Plate 80. Results from other scenarios can be seen in the Supplement file temperature.pdf,and a movie of the Base scenario can be seen at http://geosci.uchicago.edu/∼ archer/ spongebob active/temperature.active.movie.gif.
Fig. 4. (a) Total sediment accumulation rate and (b) redeposition fraction for the Base scenario.

Fig. 5 .
Fig.5.POC results from the Base case, with 12 cases as labeled plotted in Fig.5Supplement.Contours are respiration rates.Results from other scenarios can be seen in the Supplement file poc.pdf, and a movie of the Base and Bumpy simulations can be seen at http://geosci.uchicago.edu/∼ archer/spongebob active/poc.active.movie.gif.

Fig. 6 .
Fig. 6.(a) Sections of Darcy vertical velocities (flow relative to the grains), and (b) velocities relative to the sea floor, and (c) plots of w Darcy and w seafloor at the sea floor, for the Base scenario.Snapshots of other cases are shown in the Supplement files w darcy.pdfand w seafloor.pdf,and values at the sediment surface compared in w.pdf, and movies of w seafloor from the Base and Bumpy scenarios can be seen at http://geosci.uchicago.edu/∼ archer/ spongebob active/w darcy.active.movie.gifand velocities relative to the sea floor at http://geosci.uchicago.edu/∼ archer/spongebob active/w seafloor.active.movie.gif.

Fig. 11 .Fig. 12 .
Fig. 11.Methane hydrate concentration, percent pore volume, for the Base scenario with various other scenarios in the Supplement file hydrate.pdf.Solid black line is the stability boundary.Plotted with various other scenarios in Fig. 11 (Supplement).Animations can be seen at http://geosci.uchicago.edu/∼ archer/spongebob active/hydrate.active.movie.gif.

Fig. 13 .Fig. 14 .
Fig. 13.Carbon isotopic compositions compared with measurements from Pohlman et al. (2009).Short dashes are near the toe, long dashes intermediate, and solid lines are closest inshore.Heavy lines are data; thin lines are model results.

Fig. 16 .
Fig. 16.(a) Time-dependent evolution of the scenarios showing variability due mostly to the coarse grid resolution.(b) Inventory of methane hydrate at the ends of the simulations for all model scenarios.

Table 1 .
Summary of the configurations of model scenarios.Variation in the defined labile fraction of POC, 10 % and 100 % respectively, within the context of the Base case which takes 50 %.Effect of ocean temperature, changes of −2 • C, 2 • C, and 4 • C. Bubble redissolution scale height of 100 m and 2 km, relative to the Base case of 500 m.Plate 100, 80, 20 Plate subduction velocity of 100, 80, or 20 mm yr −1 instead of the Base case of 40 mm yr −1 .The slowest plate simulation was spun up for 20 Myr instead of Base 10 Myr.