References

Abstract. We present a new geologic-time and basin-spatial scale model of the continental margin methane cycle. The model, SpongeBOB, is used to simulate evolution of the carbon cycle in a passive sedimentary continental margin in response to changing oceanographic and geologic forcing over a time scale of 200 million years. The geochemistry of the sediment column is altered by the addition of vertical high-permeability channels intended to mimic the effects of heterogeneity in the real sediment column due to faults, and produces results consistent with measured pore-water tracers SO42− and 129I. Pore water dissolved inorganic carbon (DIC) concentrations are consistent with chemical weathering (CaCO3 formation from igneous rocks) at depth within the sediment column. The carbon isotopic composition of the DIC is consistent with a methane production efficiency from particulate organic carbon (POC) of 50%, which is somewhat lower than redox balance with the H / C of organic matter in the model. The hydrate inventory in the model is somewhat less sensitive to temperature than our previous results with a one-dimensional model, quite sensitive to reasonable changes in POC, and extremely sensitive to the ability of methane bubbles to rise within the sediment column, and how far gas-phase methane can get through the sediment column before it redissolves when it reaches undersaturated conditions. Hydrate formation is also sensitive to deep respiration of migrating petroleum. Other phenomena which we simulated had only a small impact on the hydrate inventory, including thermogenic methane production and production/decomposition of dissolved organic carbon.


Introduction
Models of methane hydrate formation in sediments of the deep-sea have tended to follow the example of early diagenesis modeling in adopting a one-dimensional formulation (Burdige, 2011;Chatterjee et al., 2011;Davie andBuffett, 2001, 2003a, b;Garg et al., 2008;He et al., 2011;Liu and Flemings, 2007;Xu, 2004).For the top meter of the sediment column (Archer et al., 2002), the 1-D approximation is a good one because vertical diffusion is far faster than any possible lateral pore fluid transport or diffusional effects.The methane cycle challenges a 1-D formulation, however, because the relevant chemical processes occur hundreds of meters below the sea floor, and clearly show the impact of pore fluid and gas phase migration in more than one dimension.One-dimensional models of the upper sediment require fluid flow imposed as a bottom boundary condition, and the hydrate inventory of the model is very sensitive to this parameter (Buffett and Archer, 2004).The methane cycle is also impacted by time-dependence of the sediment and organic carbon deposition rates, which are governed by onshoreoffshore processes, another benefit of incorporating the second dimension into the model formulation.
SpongeBOB is formulated in the onshore-offshore dimension laterally, and to bedrock in the vertical to internalize the entire coastal margin carbon and methane cycles, of which the methane hydrates near the sea floor are just one manifestation.Many of the processes and mechanisms, such as all the different means of downslope sediment transport or pathways for subsurface fluid flow, are impossible to represent accurately given our present state of knowledge; however, by representing them as best we are able, we hope to gain some Fig. 1.SpongeBOB in isostatic equilibrium.Colors are temperatures on the computational grid of the model, traced horizontally by the black lines.The thickness of the underlying crust grades from thin ocean crust on the right side of the domain into a thicker (and less dense) continental crust on the left.The transition between ocean and continental crust is smoothed over a spatial scale of 25 km.Local isostasy applies, on a time scale of 10 000 yr and with spatial smoothing over 20 km.An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/fig1.atlantic.movie.gif.qualitative picture of the relative importances of the various driving factors and processes.
The model is set up to simulate an accumulating passive continental margin, such as the east coast of the United States, which has been accumulating sediment since rifting of the Atlantic about 200 million yr ago.One of the beststudied locations for methane hydrates is the Blake Ridge, a "drift deposit" ridge jutting out into the Atlantic (Borowski, 2004).The Blake Ridge owes its existence to the hydrodynamics of particle sedimentation, presumably driven by existing topography and the energetic western boundary flow of the Atlantic.Although Blake Ridge is attached to the continental slope of the US, our model results (which only represent the slope, with no drift-deposit ridge) suggest that Blake Ridge might not be entirely typical of slope sediments in general, and thus might be deceptive basis for global extrapolation.
We gauge the sensitivities of this complex model by showing a variety of sensitivity model runs, varying critical parameters or turning processes on and off.The characteristics of the runs will be described as they come up in the model formulation section, and they are summarized in Tables 1 and  2.

Computational grid
SpongeBOB is formulated on a free-form grid sigma-type vertical coordinate system in depth (Fig. 1), in which the grid cells expand to fill the sediment column as it expands through time, stretching from bedrock to the sediment surface.As sediment accumulates at the sea surface each time step, sediment is also moved between subsurface grid cells in order to maintain the stretched grid.At the beginning of the simulation, all grid cells contain enough solid material to fill a thickness of 1 m.The initial porosity of the cells is "drained", such that the sediment load is carried entirely by the solid matrix with no excess pore fluid pressure to drive a Darcy flow.Sediment is deposited on the top face of the top cell, and a material flow rate is calculated for the bottom face of the top cell, such that the top cell keeps some of the new material for itself but passes most of it on fill the subsurface cells below.The process continues through the sediment column, and the boxes all expand uniformly.
A few grid cells in the top of the domain are restricted, such that they cannot expand beyond a thickness of 200 m, in order to maintain a minimum computational resolution in the region of the hydrate stability zone.After enough model time has passed that a restricted cell reaches its maximum thickness, it passes along the entire incoming advective flux from above into the cell below.

Isostasy
The sediment column elevation is governed by isostasy, the idea that the crust and all the material sitting on top of it are all floating in a denser fluid of mantle material (Watts, 2001).Laterally, the model domain spans between continental crust and ocean crust, as can be seen in the transition in the thickness of the crust in Fig. 1.
The effect of cooling with age on the elevation of the ocean crust is driven by the evolution of a time-dependent thermal boundary layer, extending through the crust into the mantle below.As the material in the crust and the mantle boundary layer cool, they contract and become denser.The thickness of the boundary layer is assumed to increase according to z boundary = √ κT (Parson and Sclater, 1977) where κ is a thermal diffusivity.The temperature profile with depth is assumed to be linear allowing us find the temperature at the crust/mantle boundary by interpolation, from which to compute the average temperatures of the crust and the mantle boundary layer.The temperature at the base of the mantle boundary layer is taken to be 1600 K (Iamori et al., 1995).The floating pile in the ocean crust end-member consists of the sediment column, the ocean crust, and the thermal boundary layer as it thickens and thermally contracts through time.
The continental crust is simpler because there is no timedependent cooling.The crust is defined to be 13 km thick and

No Thermogen
Thermal methanogenesis disabled Decreases hydrates by about 50 % its density is 3.0 g cm −3 .Any sediment load also contributes to the mass of the column, and the equilibrium elevation of the bedrock.The blending of domains is done by calculating pure endmember continental and ocean crust scenarios at each model grid point, and then scaling between them using an ad-hoc ocean fraction variable, ranging from 0 for continent to 1 for the ocean end member.The transition is set up with a horizontal length scale of 50 km.For each end member, continental and oceanic, the masses and heights of the column are computed.The two scenarios are blended together at this stage, taking a composite column mass that scales between the two end members, and a composite column height.Although this is not a mechanistic simulation of crustal dynamics at a passive margin, it is sufficient to produce a final envelope of sediment that resembles seismic sections of passive margins (Fig. 4), as a driver for the sedimentary core of the model.
Sea level is defined as 4000 m elevation above generic ocean crust.The hypothetical elevation of the top of the mantle, the sea in which the crustal material floats, is defined to be consistent with this.
The equilibrium heights for the top of bedrock are smoothed horizontally through the grid using a diffusion algorithm that achieves a spatial scale of a few tens of kilometers (three passes with a diffusion coefficient of 10 6 m 2 yr −1 ).Flexure and rigidity of the plates becomes important in accretionary margins, but for this study the lateral interaction in isostasy is rudimentary, merely setting the stage for the sediment column.The elevation of the bedrock evolves toward the isostatic equilibrium value with a time constant of 10 −4 yr −1 .

Deposition
The goals of the sedimentation scheme are to create appropriate sedimentation rates and grain sizes, in adaptive response to a time-evolving submarine landscape.The existence of a shelf break, and the prograding clinoform pattern of sediment deposition, require the existence of a region of fast sediment accumulation just off the shelf break, a sediment depocenter (Pirmez et al., 1998).Sediment input from the continent is also sorted by grain size, with largest grains depositing closer to shore and finer grain sediments offshore (Fig. 2).
Offshore sediment transport and deposition in Sponge-BOB is parameterized in an approach following Pirmez et al. (1998).Sediment comes into the model domain on the continental (left-hand) boundary, suspended in the ocean.It advects offshore with a sediment transport velocity, u ocean , that decreases inversely with water depth, as if the ocean margin were a giant river carrying suspended particles away from the continent in a plug flow.Suspended particles settle out following Stokes' law according to their grain size where R Particle is the particle radius, ρ stands for density, and g is Earth's gravitation acceleration.Particles deposit only if the shear stress at the bed, calculated as is too low to lift the particles faster than they sink, determined as a critical shear cutoff The sediment transport scheme is not mechanistically realistic, but is intended to serve as a realistic boundary condition driver for the sediment column core of the model.The sediment transport offshore velocity (0.1 m s −1 in water depth depth of 10 m) and the resulting shear deposition scheme were tuned to produce a reasonable facsimile of a seismic section of the Atlantic margin (Fig. 3).
The critical shear velocity depends on particle radius, as does the settling speed, effectively sorting the particles by size with offshore distance.The impact of water depth on the plug flow and hence the shear velocity allows changes in sea level to impact sedimentation and the evolution of the sediment column.

Slope-driven offshore transport
Turbidites in abyssal sediments attest to the importance of slope failure to the margin sediment dynamics (Meiburg and Kneller, 2010).These processes are handled as a function of the sea floor grade.If the grade exceeds a critical value, here taken to be 6 %, the material that would deposit according to the shear/sinking speed criterion is instead classified as "resuspended".Also, if the sea floor slope ever manages to exceed the critical slope, surface sediment is eroded into the resuspended pool.The resuspended material redeposits using the same settling speed/offshore transport criterion as used for the primary settling material, but resuspended material only begins to deposit where the sea floor slope is less than 1 %, reflecting the difference between a turbidity flow which only stops on the plain as opposed to individual sinking particles assumed in the primary sedimenting material.The particulate organic carbon (POC) concentration of the sediment is determined in the model as a function of water depth at the time of initial sediment deposition, as discussed under the heading Carbon Cycle.When material is resuspended from the sea floor, in contrast, its POC concentration is conserved and redeposited.Because the offshore transport is assumed to be turbidites, we assume it escapes the degradation that individual settling particles, the primary deposition, are subject to as they slowly sink through the water column.Resuspended organic carbon is assumed to be distributed among the grain size fractions of the resuspended material according to surface area, selectively concentrating (according to mass) the organic carbon into the smaller size fractions.
The sediment transport scheme is intended to be a utilitarian parameterization of sediment transport processes, rather than a mechanism-resolving model such as Sedflux (Syvitski, 2001).The scheme as we have tuned it achieves the goal of a well-defined shelf break, driven by a strong sediment depocenter just off the shelf break, and response to sea level changes that mimic the Sedpak model (Strobel et al., 1989).The accumulation rate at 2 km depth in the model is close to that of Blake Ridge, about 2 cm kyr −1 .Total sediment thicknesses and extents are comparable with a seismic image of the Atlantic margin from Kennett (1982) (Fig. 3).Model sensitivity to these sediment transport processes was gauged by altering the critical sea floor slope in sensitivity simulation Lo Slope.

Sea level changes
Sedimentation is heavily impacted by eustatic changes in sea level, with highest rates close to shore but in water that is not too shallow.Sea level changes are imposed on the simulation as a 240 million yr cycle, peaking at high sea level about 100 million yr ago and declining since then, to generally and smoothly follow the larger-scale sea level changes in the Haq et al. (1987) curve (Fig. 4).The organic carbon (POC) of the  new sediment covaries with sea level in the model (see Carbon Cycle section).Two sets of simulations were done, one which neglected time-dependent sea level changes (Table 1), and another set which included these effects (Table 2).Results in this paper will be presented as time slices from various points in the simulations, but animations of the entire runs are available as links from the figure captions.

Porosity and excess pressure
The pore fluid excess pressure, which drives the pore fluid flow, is calculated from a constitutive equation in which the excess pore fluid pressure can be calculated from the porosity and the mass load overhead (Flemings et al., 2002).The idea is that there is a one-to-one relationship between the porosity of the sediment and the amount of solid load it can support.The compression of the bulk sediment by exclusion of pore water at hydrostatic pressure is governed by where the β-values are bulk compressibilities and S supp is the solid load that can be supported by sediments of that porosity.Since pore fluids are relatively close to hydrostatic, the porosity profile in the sediment column is largely governed by this relationship.We find that to produce the relatively high porosity of surface sediments and the lower porosities of the deep sediment, the column requires two exponentials with different β values (Fig. 5b).Coefficient β 1 is 0.0312 and produces compaction on a depth scale of about 10 km, representing ductile flows and recrystallization of mineral grains.The value of β 2 is 0.2, producing a compaction depth scale of about 2 km (Fig. 5), representing the collapse of clays in the shallower sediment column.Because of the double exponential form of the constitutive equation, calculating S supp from is an iterative calculation.The excess pressure is then calculated as where S v is the total load of the sediment column from all phases, and P hydro is the hydrostatic pressure (the load due to the fluid phase).

Temperature and heat flow
Subsurface temperatures are primarily governed by thermal conduction, a diffusional process, although fluid advection and latent heat also affect the temperature.The thermal conductivity and heat capacity are taken as a volume-weighted averages of the components: sediment, pore fluid, methane hydrate, and methane gas.A geothermal heat flux is imposed at the base of the sediment column of 70 m W m −2 , with sensitivity runs spanning from 50 to 120 m W m −2 .

Clay dewatering
At elevated temperatures, hydrated clays such as montmorillonite can release water to form illite. Montmorillonite is taken to be the dominant sedimenting material , and it dehydrates at temperatures above about 50 • C (Yang, 2006) with an activation energy of 196 KJ mol −1 .The total volume of water released in the model by this mechanism is probably an overestimate, and because the rate constant is poorly tuned in the model, the dehydrate takes place at a depth which is shallower (colder) than reality.The impact of this process is gauged in a simulation in which it is neglected called No Clay Dewatering.

Permeability
The vertical permeability k v is parameterized as where r is the mean grain radius and φ the porosity (Bear, 1972).A factor of 100 anisotropy is imposed, following the example of Daigle and Dugan (2011), such that to parameterize for the effect of layers like sandy turbidites, that are unresolved in the model.We show results of two model permeability formulations, one as presented so far, and the other with added vertical channels of enhanced permeability, such as faults, pipes, channels, and similar structures (cases Base and No Channels).The channels are situated every 5 grid cells throughout the model domain, and remain stationary as the continental margin progrades offshore through the simulation.The horizontal grid spacing is 3.15 km, and the vertical permeability in the channel is 10 times higher than given above.This is only a crude representation of a realistic flow regime through faults and permeable channels that seem to dominate subsurface flows (Flemings et al., 2003).Most of these features would be much smaller than the grid resolution of the model can represent more realistically, although some seismic "blowout zones", disturbance of the layered sedimentary structure due presumably to gas escape from the deep sediment column, do span several kilometers in horizontal extent (Hill et al., 2004).The intent is to gauge in a qualitative way the sense and potential magnitude of the impact of heterogeneity on the flow and chemical evolution of the sediment column.
The permeabilities and the excess pressures of the two simulations at the end of the simulations are shown in Fig. 6.

Darcy flow
Fluid flow follows Darcy's law in response to the pore fluid excess pressure gradient (Athy, 1930), to find along-and across-grid velocities as where µ is the dynamic viscosity, k the permeability, P excess the excess pressure, and x and z the grid dimensions.All of these variables except the viscosity are defined on the two grid dimensions i and j .The velocities are defined on the cell faces.
A flow limiter in the vertical protects the calculation from numerically overstepping (for thin, permeable layers) by testing for the change in fluid pressure that would result from unimpeded Darcy flow each time step: where the derivative of the porosity/S supp relation ship (Eq. 1) is stored for each grid cell as part of the iterative solution of Eq. ( 1) for S supp .If the vertical pressure gradient would change by more than 10 % in a time step, a stand-in vertical velocity is calculated that limits the change to 10 %: The idea is that a thin sandy layer would quickly drain, reaching zero excess pressure gradient, and the stand-in velocity is designed to achieve that end in a stable way even if the time step for Darcy flow would need to be much shorter.The 10 % coefficient implies that it will take approximately 10 time steps (5 yr) for the layer to relax to a drained condition (any faster leads to numerical instability).Most of the vertical flow in the model however is well-resolved in time, leaving the flow limiter for a few "emergency" grid points and times.

Total fluid flow including sedimentation and grid stretching
Darcy flow is defined relative to the solid grains.This is a different definition of fluid flow from the convenient and usual metric in sediment diagenesis models, which is defined relative to the sediment-water interface, which we will call w seafloor .Formation of methane hydrate in models is very sensitive to the value of w seafloor , upward carrying methane into the stability zone from below, or downward away from the stability zone (Buffett and Archer, 2004).At the sediment surface, the two definitions of vertical velocity differ because the sediment surface is moving as solid and fluid material accumulates.The burial of fluid in the pores between the sedimenting solid constitutes a vertical flow across the sea floor, when the reference frame is the moving sea floor itself.At the sedient surface, where SedRate is the bulk sedimentation rate and φ 0 the porosity of freshly sedimented material.
Beneath the sea floor, the calculation of w seafloor requires two additional terms.One is the movement from one grid where ∂z boundary /∂t represents the movement of the cell boundary.
For subsurface grid points, w seafloor is calculated as The third term represents the relative movements of the seafloor and the depth of the subsurface grid cell in question.Like w bur and w Darcy , w seafloor is defined on the cell faces.Numerically, the velocities are combined before any fluxes are calculated, using the net flow across a grid boundary within a given time step, rather than advecting material back and forth in multiple operations which would result in needless numerical mixing.
A subsurface value of w seafloor is defined relative to the moving sea floor, so that if a subsurface parcel were flowing upward just quickly enough to remain at a constant depth below the seafloor, its w seafloor would be defined as zero.Sections of w Darcy and w seafloor are shown in Fig. 7, and values at the sediment surface in Fig. 8.
When vertical permeable channels are added to the simulation, the channels take most of the compaction flow.The horizontal anisotropy applied to the permeability allows the sediment columns in between the channels to vent fluid into the channels for expulsion.Interestingly, for the formation of methane hydrate discussion to come, the channeled simulation finds strongly downward w seafloor in the grid cells in between the channels.These cells have vented their upward Darcy flow into the permeable channels, leaving the burial flux of fluid with sedimentation to drive w seafloor in the downward direction, where it brings methane-depleted fluid into the hydrate stability zone from above, rather than methanerich fluid from below.

Deposition
The carbon cycle of the deep biosphere is driven by the rain of particulate organic carbon (POC) to the sea floor.Organic carbon deposition and concentrations in surface sediments of the ocean depend on primary productivity in the ocean, runoff of terrestrial POC in some places, degradation in the water column, offshore transport, and degradation in the surface sediments which depends on the oxygen concentration of the overlying water and perhaps the temperature.Rather than mechanistically simulating that complexity, the POC concentrations and characteristics of the sedimenting material in these simulations are determined as a simple function of water depth of initial deposition and time, allowing us to control POC deposition as a boundary condition to the model rather than simulating it as a prognostic output.We used the POC parameterization to drive a scenario in which periods of high sea level have enhanced carbon preservation, a cartoonish representation of the transition from the Cretaceous to the present (Hunt et al., 2002).
Although the link between ocean oxygenation and POC deposition is not mechanistically clear or straightforward to model, we concoct and drive the model with an ocean oxygenation "index" that varies from 1 to 4, with 1 the most oxic and 4 anoxic.Ocean state 1 is called "Atlantic", and the sedimenting POC concentration profile is comparable to sediment surface POC concentration off Cape Hatteras today (Mayer, 2002).In this state, the POC concentration is at most 2 % dry weight POC, with an H/C ratio of 0.7.The anoxic end-member buries up to 5 % POC, with an H/C ratio of 2 (Fig. 3b).
The ocean oxygenation index covaries with sea level and the deep ocean temperature, resulting in a scenario in which a warm high-sealevel high POC deposition ocean (100 million yr ago) transitions into the characteristics of the Atlantic ocean by the end of the simulation.For a given grid cell and time step, the POC concentrations are interpolated in depth, then interpolated in the ocean state, so that carbon deposition can vary smoothly in time, in synchrony with changes in eustatic sea level.
The impact of a persistent change in ocean oxic state on hydrate abundance was gauged in simulations with timeuniform forcing, called Atlantic, Pacific, OMZ, and Anoxic.In the time-dependent simulations, high POC concentrations accumulate in a time of high sea level in the simulation, which becomes a buried feature by the end of the simulation (Fig. 9 left).

Biological degradation
Biological degradation is taken to be a function of temperature (T , • C), sediment age, and dissolved methane concentration.
where k resp is a first-order rate constant for POC degradation.The activation energy Ea is 70 kJ mol −1 , resulting in a tripling in the rate with a temperature increase of 10 • C.This is a relatively moderate temperature sensitivity; Burdige (2011) argues that it could be 200 kJ mol −1 , which would increase the reaction rate by a factor of 24 for an increase in temperature by 10 • C.
The second term represents inhibition of biological activity above a critical temperature of 50 • C. The formulation uses a smoothed "Heaviside" function to shut down the respiration reactions above a threshold temperature of 50 • C.This function is used several times in the SpongeBOB formulation, defined as This function provides a smooth transition between the on and off states (1 and 0), rather than the discontinuous step Raising the function to a power, as done here, tightens the transition, here intended to mimic a fairly strong shutdown of respiration rates when it gets too hot.The third term represents the correlation between carbon age and its reactivity, mimicking the astonishingly wide correlation observed in nature (Hedges and Keil, 1995;Middelburg, 1989).The correlation between age and reactivity does not imply a direct or simple causation, and "priming experiments" (Bianchi, 2011) demonstrate that the age parameterization cannot be universally true.But the correlation between age and reactivity captures many orders of magnitude of variation, and seems essential to account for (Wallmann et al., 2006).
The effect of the "metabolite inhibition" effect as proposed by Wallmann et al. (2006) is represented in the last term in the formulation.Wallmann included dissolved inorganic carbon concentration in addition to methane, but for the sake of model stability and development we took the simplifying liberty of excluding DIC inhibition of methanogenesis.
Respiration is concentrated in the top 2 km of the sediment column, limited by temperature (Fig. 10).It is also generally focused offshore of the shelf break, except during periods of very high sea levels, during which POC can accumulate in the shelf sediments.When sea level falls, the biologically available POC in the shelf sediment burns out, leaving most of the biological activity in slope sediments, and in particular in the sediment depocenter just off the shelf break.
Model sensitivity to the respiration formulation was gauged by varying the "labile fraction" of the POC (cases Bio 10 %, Base (which is 50 %), and Bio 100 %).

Biological methanogenesis
The respiring carbon is assumed to produce DIC and molecular hydrogen, a transient reducing agent.If sulfate (SO 2− 4 ) is present in the pore water, the hydrogen reduces and consumes it; otherwise, it reacts with CO 2 (dissolved inorganic carbon or DIC) to produce methane in reactions where the variable x denotes the relative concentration of POH/POC in the reacting organic matter.
The maximum efficiency with which the biologically produced CO 2 can be converted to CH 4 is determined by conservation of carbon, hydrogen, and oxygen from the organic matter to methane plus CO 2 , realizing that supplemental oxygen and hydrogen comes from H 2 O.This maximum CH 4 production efficiency is which works out to 62 % if the H/C ratio of the reacting organic matter is 1.However, carbon isotopic systematics and data, described below, will restrict the efficiency of CH 4 production somewhat, as if about 20 % of the H 2 reacts with oxidized mineral phases, such as Fe 2 O 3 rather than with CO 2 to make methane, leaving an overall methane production efficiency of about 50 %.Methane is oxidized in the presence of SO 2− 4 by the anaerobic oxidation of methane (AOM) reaction.The reaction in nature is fast enough to permit just a very narrow depth range of overlap between the high-SO 2− 4 pore waters in the surface sediment column and high-CH 4 waters below (Burdige and Geiskes, 1983;Chatterjee et al., 2011;Davie and Buffett, 2001).The reaction in the model is based on whichever species is stoiciometrically limiting, with a effective rate constant of 5 % yr −1 , unrealistically slow but numerically stable and able to eliminate simultaneous coexistence of SO 2− 4 and CH 4 in any grid cell.SO 2− 4 is only found in a few grid cells just inside the continental shelf break.The grid resolution of the model is too coarse to resolve the details of pore water SO 2− 4 , so we leave an analysis of the behavior and impacts of SO 2− 4 to future work.

Petroleum and DOC
The model produces petroleum based on chemical and thermal conditions.Petroleum formation in the real world is limited by the H/C ratio of the organic matter, proceeding if the H/C is greater than 1, and continuing until the ratio approaches 1.The temperature dependence of the rate constant is governed by an activation energy of 150 kJ mol −1 , with a pre-exponential constant of 1 × 10 16 yr −1 .These parameters confine petroleum generation to a temperature window of 60-150 • C, as inferred from the distribution of oil (Hunt, 1995).Petroleum genesis temperatures are too warm to allow bacterial consumption, but we gauge the potential impact of bacterial respiration of petroleum by assuming that 10 % of the oil is upwardly mobile at a rate of 1 m per thousand years.If it reaches cool enough conditions it is consumed by bacterial respiration, following the temperature and solute concentration dependence (but not the POC age dependence) of bacterial POC degradation kinetics above.The rate of methane production by this process of "secondary respiration" is shown in Fig. 10b to be deeper than POC respiration, which is more confined to the upper sediment column.The model also assumes "sloppy feeding", in which 10 % of biological POC consumption goes to dissolved organic carbon (DOC).Simulations to determine the importance of these processes in the model are called No Petro and No DOC.

Thermal methanogenesis
Methanogenesis from thermal degradation of organic carbon proceeds according to an Arrhenius equation with an activation energy of 230 kJ mol −1 and a pre-exponential rate constant of 4 × 10 21 yr −1 , which produces a lifetime of a million years at about 160 • C. In contrast to the biological rate of methane production, the rate constant for thermal methane generation is very well constrained and consistent in field data, driven by the hunt for petroleum (Hunt, 1995).
The rate of methane production is taken to be this rate constant multiplied by the concentration of particulate organic hydrogen, POH.The POC inventory decreases by 1/4 of the decrease in POH during methanogenesis.The distribution of  and D).An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/ch4.movie.gif.
the H/C ratio in the sediment column (Fig. 9 right) is primarily driven by the effects of changing ocean chemistry (anoxia) and thermogenic degradation.The H/C ratio approaches zero as the carbon matures, and the decrease in H directly limits the production of methane in the model.

Methane bubbles and hydrate
The dissolved methane concentration is compared with the solubility of dissolved methane, relative to the two possible other phases of gas bubbles or hydrate, in Fig. 11.Hydrate or bubbles both release or take up methane from the dissolved phase.Disequilibria in the three-phase methane system is restored toward equilibrium by conversion to or from bubbles or hydrate with a numerically-tractable but ad-hoc time constant of 1000 yr.This unrealistically slow time constant, a compromise that should doubtlessly be higher for future detailed simulations, permits the dissolve methane concentration to reach unrealistically high supersaturation values ( = 1.2 in Fig. 11), but since dissolved methane is a transient intermediate species in the pathway POC to CH 4 to bubbles, the rate of bubble formation will be determined by methaneogenesis, and the dissolved methane concentration in the steady state will grow in to fulfill that balance, leaving only the question of whether the altered CH 4 concentration changes some other aspect of its dynamics, by diffusing more quickly for example.
When the precise depth of the stability boundary (determined by depth interpolation) falls within a grid cell, the nongaseous CH 4 is allocated into bubbles and hydrate in proportion to the volumes of the two zones, and the concentrations of the phases, for purposes of advection and diffusion, are calculated using the partial grid cell volumes.This prevents  E and F).An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/bubbles.movie.gif.
a numerical artifact that otherwise appears as the depth of the stability boundary crosses through the coarse numerical grid of the model.The methane cycle in the SpongeBOB requires gas-phase mobility in the sediment column.If bubbles are not allowed to migrate, gas completely fills some parts of the sediment column pore space (Fig. 12c).Gas volume in the pore space of the model currently and unrealistically has no impact on the fluid pressure or flow, allowing unlimited gas accumulation to fill up to 200 % of the original pore volume.That much methane simply has to go somewhere.
In a basin-scale sense, gas mobility comes as no surprise.The view from petroleum geology is that, on geologic time scales, sediment columns are permeable to gas except under permafrost soils and evaporites, which overlay the largest gas accumulations in the world in Siberia and the Middle East, respectively.Other parts of the world, for example the North Sea, are vastly enriched in oil relative to the expected amount of gas were it retained (Hunt et al., 2002).On the smaller spatial scales of structural hydrate deposits, such as at Hydrate Ridge (Treguer et al., 1995), gas migration is also apparent.Gas migration is a fundamental part of the Liu and Flemings (2007) model to explain the seismic "wipeout zones" and observed methane emissions from the seafloor within the hydrate stability zone (Hill et al., 2004).
In contrast, for modeling "stratigraphic" methane hydrate deposits typical of most of the continental margin, gas migration is typically neglected or assumed to be small (Burdige, 2011;Chatterjee et al., 2011;Davie and Buffett, 2001;Davie and Buffett, 2003a;Davie and Buffett, 2003b;Garg et al., 2008;He et al., 2011;Liu and Flemings, 2007;Xu, 2004).These one-dimensional methane hydrate models span a domain a few hundred meters below the stability zone boundary and assume that the bubbles there are immobile.This is not a priori wrong; it could be that gas transport is only important in the structural hydrate deposits, more obviously associated with faults and channels.However, we do find that methane hydrate formation in SpongeBOB is strongly driven by upward migration of gas, the same gas flux as is required to evacuate excess gas from the deep sediment column.
Gas tends to form if the dissolved concentration of the gas exceeds the Henry's law solubility, a function of the pressure of the gas.Surface tension raises the pressure of the gas phase somewhat over that of the pore fluid, but the differences are small relative to the overall pressure, so the impact of excess gas pressures on methane solubility have been neglected in the model.Bubbles typically grow, in cohesive sediments, by fracturing the sediment matrix (Boudreau et al., 2005;Jain et al., 1995), requiring further excess pressure that is a function of sediment properties such as fracture strength.With a continuous source of methane, dissolved methane concentrations can rise until Henry's law produces gas at whatever pressure is required to create the bubble.
It is another issue to understand the upward mobility of the gas, as opposed to its growth in place.Even at pressures near the base of the sediment column, the compressed gas is buoyant, manifesting as an upward body force that intensifies as the bubble rises and decompresses.One source of excess pressure at the top of a bubble is the hydrostatic pressure difference between the top and the bottom of the bubble (Jain and Juanes, 2009).However, for a bubble several kilometers deep in the sediment column to break through the sediment column by this mechanism, the bubble would have to several kilometers tall, to get the excess pressure at the top of the bubble high enough to break through the sediment column.Presumably gas migration from the deep sediment resembles gas blasting through a pipe more than it does a long bubble contemplating rising.The gas can manufacture its own means of escape by building up to lithostatic pressure, possibly venting episodically.
Laboratory experiments find a critical gas saturation (pore volume fraction) of about 10 % where gas begins to mobilize (Leas et al., 1950), and values of 1-2 % are often assumed in models for natural gas in wells (Reagan and Moridis, 2007).Higher gas saturations connect bubbles together, crossing the hydrostatic pressure gradient as well as providing a pathway for the gas to flow, which is important because gas is much less viscous than pore water, even when it is compressed.The standing stock of bubbles in the sediment column is another transient intermediate reservoir for methane in the pathway POC → dissolved methane → standing gas bubbles → migrating gas.In the steady state, the standing gas bubble reservoir will find the concentration at which the methane fluxes balance.The essential balance is between methanogenesis and upward gas flux in the model.A bias in the bubble flux parameterization would result in an incorrect standing bubble concentration but no error in the bubble flux itself, unless the bubble concentration leads to some distortion of the methane budget, such as by burial of gas.
We have adopted a parameterization of bubble migration that depends on the bubble fraction in pores CH 4 loss=H smooth (F b −10 %, 5 %) 3   where again the H smooth is a smoothed version of a Heaviside function.The bubbles escape as their concentration approaches or exceeds a critical pore volume fraction of 10 %.In practical terms, the bubbles in SpongeBOB only fill at most 1-4 % of the pore space, because the loss rate, however slow, is able to keep up with production, even when the saturation value is much less than the parameterized critical value.
Rather than resolve the rise of the gas through time, we instantaneously redistribute some of the CH 4 lost as migrating bubbles into the dissolved phase to grid cells above where the bubble first started moving.This approach is analogous to the treatment of sinking plankton in the water column of many ocean chemistry models, in which the sinking organic matter is assumed to reach some mean distance before instantaneously redissolving.The bubbles once in motion rise through the sediment until they encounter pore water that is undersaturated with respect to dissolved methane, at which point the upward flow is attenuated with height as CH 4 flux(z)=CH 4 flux(0) • e −z/z scale where the value of the scale height parameter z scale is taken to be 500 meters in the high-resolution run, but varied to 100 and 2000 m in two of the sensitivity runs (called Bubb 100m and Bubb 2k).The methane cycle in the model in general, except for the standing bubble inventory, is not very sensitive to the bubble mobility parameterization, but the methane hydrate inventory is extremely sensitive to what happens to the gas after it begins moving.
The distribution of bubbles resulting from these assumptions is shown in Fig. 9.A standing crop of bubbles is sustained by respiration in the upper sediment column, and there is a column of bubbles carrying methane from the thermogenic zone up toward the sediment surface.Thermogenic methane in the model systematically evades the hydrate stability zone, because the sediment column is only deep enough to reach high enough temperatures when the top sediment column has nearly reached the sea surface, so that there  (B) Base case profiles compared with data in red from Blake Ridge ODP sites 994 and 997 (Rodgiruez et al., 2000)and the West African margin Sivan et al. (2007), and (C) with Urey reactions that consume CO 2 disabled (note change in color scale from A).An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/dic.movie.gif. is in general no hydrate stability zone over thermal methane production in the model (Fig. 9).This geophysical mechanism could lead to a general bias toward biogenic rather than thermogenic methane in deep margin hydrate deposits.

DIC and subsurface weathering
Measurements of the titration alkalinity of pore waters (a reasonable facsimile to DIC under pH-neutral conditions) from Blake Ridge reach fifty times higher than the ocean value within the top few hundred meters, but then either decrease to half that by the deepest samples measured, about 700 mbsf (Rodgiruez et al., 2000).A similar decrease with depth was documented from the West African margin (Sivan et al., 2007).
In the model, if CO 2 is allowed to accumulate in deep sediment pore water without a compensating sink, the DIC concentration of the model increases monotonically to bedrock (Fig. 13b and c).The difference can be explained as a chemical weathering or carbonate formation reaction (Wallmann et al., 2008), such as anorthite + CO 2 → kaolinite + CaCO 3 .
Throughout the upper sediment column, the equilibrium concentration of dissolved CO 2 based on this equilibrium is much lower than the measured concentrations.At low temperatures typical of the surface of the Earth or the upper sediment column, this reaction is too slow to reach

D. E. Archer et al.: A model of passive margin carbon and methane
equilibrium, but at a high enough temperature the reaction proceeds quickly as the metamorphic decarbonation reaction.We treated this reaction as a first-order uptake in DIC.This neglects the impact of the alkalinity or pH of the solution on the CO 2 concentration, but the constancy of porewater alkalinity/DIC ratios (Sivan et al., 2007) and pH make this seem a reasonable assumption.The reaction as written will affect the concentration of calcium, Ca 2+ , but it is not calcium limited unless there is a shortage of unweathered igneous calciumbearing rock.The model assumes equilibrium with CaCO 3 , rather than try to resolve the porewater chemistry into its components alkalinity, calcium, and DIC.
The kinetics of DIC uptake are taken to be a function of temperature as where the activation energy Ea is 120 kJ mol −1 , increasing the rate constant by a factor of five in response to a temperature increase of 10 • C. The slowest pre-exponential multiplier we tried was 3 × 10 13 yr −1 , which produces a reaction timescale of about 10 8 yr at 15 • C, so that essentially nothing happens in surface sediments, to a time constant of 1000 yr at 100 • C. We also show profiles of model DIC concentration with faster uptake kinetics, by factors of 3 spanning two orders of magnitude.
The model profiles (Fig. 13b) are compared with alkalinity data for ODP sites 994 and 997 on the Blake Ridge (Rodger, 2000).Both the peak concentration and the depth of the peak are best simulated by the model with an intermediate rate constant of the range we tried, 3 × 10 14 yr −1 .
For all runs with nonzero Urey reaction rate constant, there is a drawdown in DIC concentration with depth.The similarity to measured concentrations lends strong support to the significance of this reaction to sediment pore water DIC.It is no surprise that carbonate formation occurs within the sediment column (Rodger, 2000), but it is interesting to note the impact that this process has on measured DIC concentrations.
Also in all weathering runs, the concentration of DIC in sediments warmer than 100 • C in the model is effectively controlled by equilibrium, determined from thermodynamic data for the phases in the weathering reaction.The highest DIC levels in the simulation are found in the hottest, deepest sediments due to the effect of temperature on the equilibrium condition, the beginnings of metamorphic decarbonation (Fig. 13a).At temperatures between about 50 to 100 • C, the equilibrium DIC concentration is relatively low, but the reaction kinetics in the model are still fast enough for equilibrium to be reached, pulling DIC concentrations down in that depth range, and isolating the high DIC concentrations in the hottest sediments from reaching the sediment surface.

Carbon isotopes
The model tracks parallel tracers representing the stable carbon isotopes in DIC and CH 4 .Using methodology used in ocean (Maier-Reimer et al., 1982) and sediment modeling (Jorgensen, 1978), the isotopologs undergo the same advection, diffusion, and chemical reaction operators, but with offsets in the reaction rates corresponding to the fractionation factors.Rather than try to numerically resolve the large relative concentration differences between the isotopes, the isotopic standard ratio for all isotopic systems is defined to be 1 meaning equal concentration of both isotopes at 0 ‰.The results, presented in ‰ deviations from the standard ratio, are directly comparable with isotopic ratios of natural samples presented relative to their separately defined standard ratios (PDB for δ 13 C, for example).
The isotopic boundary conditions are that the DIC in the ocean is assumed to be 0 ‰, and organic matter is assumed −25 ‰.Organic matter is respired to DIC without fractionation.CO 2 reduction carries a fractionation of −90 ‰ at 0 • C, dropping to about −55 ‰ at 50 • C (Whiticar and Faber, 1986).Thermogenic methane is assumed to have the isotopic composition of organic matter.Diffusion and anaerobic oxidation of methane (AOM) also fractionate the carbon isotopes.
The isotopic composition of DIC in the sediment column (Fig. 14) is lightened by degradation of organic matter, but then pulled heavier by extraction of isotopically light CH 4 .As a result, the isotopic composition of DIC or alkalinity in porewater samples starts near 0 ‰ at the sediment surface, then trends light in the near subsurface due to respiration in the non-methanogenesis zone, then gets heavier, reaching a maximum value of about +10 ‰ (Fig. 14).These features are mimicked by the model, but the degrading effect of the coarse grid resolution can clearly be seen in the broadening and attenuation of the subsurface negative isotopic excursion.Sivan et al. (2007) interpret the maximum δ 13 C of the DIC as an indicator of the rate of methanogenesis, but it turns out to be more sensitive in the model to the relative partitioning of respired carbon between DIC and methane.If the respired carbon has an isotopic composition of −25 ‰, and the fractionation between CH 4 and DIC at some temperature were 60 ‰, for example, then a balance of about 50 % methane versus DIC would tend to produce DIC at an isotopic composition of +10 ‰.In time, this source composition dominates the pore water value.Model experiments did not show a sensitivity of δ 13 C in DIC to the overall methanogenesis rates, but it does respond to the relative CH 4 production efficiency.The best fit seems to require some loss of H 2 in the pathway, as if by reduction of mineral phases.

Iodine-129
Iodine-129 has been used as a tracer for the age of the degrading organic matter in a sediment column (Fehn et al., 2000).Iodine-129 is produced by cosmic ray spallation in the atmosphere and it has a 22.8 million yr e-folding lifetime.Radioactive and stable isotopes of iodine are incorporated into biological material, and released to the pore water as the carbon degrades biologically.Iodine released by older organic carbon will have less iodine-129 because of radioactive decay.In the model, both stable and unstable iodine are carried with the pore flow and are subject to diffusion along with the rest of the dissolved tracers.Fehn et al. (2000) measured I 129 ages 55 million yr old in pore waters of solid sediment that was 1.8 to 6 million yr old, a difference of about 50 million yr, which they interpreted as indicative of 50 million yr old carbon degrading.
Iodine-129 seems like a very rich tracer in the Sponge-BOB model, indicative of the history of the accumulation of the sediment column and the expulsion of fluid from it (Fig. 15).Variations in I 129 age arise from compaction-driven fluid flow.When fluid flow is prevented in the model, the iodine ages correspond to the sediment ages (Fig. 15d).When fluid flow is enabled but without the permeable channels (Fig. 15b), there are iodine ages 60 million yr older than the solid ages at its maximum, located in mid-column sediments under the deep continental slope, along with I 129 ages less than the solid age in its minimum under the continental shelf.This pattern persists when permeable vertical conduits are available, with a superimposed pattern that the I 129 ages within the conduits are older than in the surrounding lowflow columns (Fig. 15c).The old I 129 ages "grow in" as the continental margin ages (compare Figures 15a and c).
The iodine ages in the depth range of Blake Ridge in SpongeBOB never get as old as Fehn et al. (2000) measured at Blake Ridge.The values in the conduits reach 20 million yr old, while Fehn measured 50 million yr.The model results  suggest that the old I 129 ages measured at Blake Ridge probably have more to do with pore fluid flow than with the age of respiring organic carbon.They also suggest that Blake Ridge pore water geochemistry may be affected by focused subsurface flow, such as crudely approximated by the model, and thus would over-represent the abundance of methane hydrate along margin sediments overall (which cannot all be the subject of focused upward flow).

Scenario runs
In addition to the "Base HR" and "No Conduits HR" cases, a series of sensitivity simulations were done, at lower spatial resolution of only 10 grid points horizontally instead of 40.
The setup, intent, and summary of the conclusions from the sensitivity runs are described in Tables 1 and 2, and the resulting methane hydrate inventories shown in Fig. 16a and b.One suite of simulations was forced by time-invariant forcing, comprising sea level, ocean "oxic" state (POC deposition patterns), and temperature with variations as described in Table 1.These allow a clearer diagnosis of the model sensitivity to these variables.Another suite of simulations imposes time-varying forcing shown in Fig. 3, as described in Table 2.
The time evolution of sea level, POC inventory, methanogenesis rates, and hydrate inventories are shown in Fig. 16c for the two representative cases, the Base case, which including the sea level changes, and case T 0, which does not.

Sensitivity to temperature
The effect of ocean temperature on model hydrate inventory is shown in Fig. 17a.These are from time-invariant suite of runs, with temperature a function of water depth only.The inventory of hydrate decreases roughly linearly with temperature, finally reaching zero hydrate with a warming of 8 • C. The active margin simulations (presented in the companion paper) have a similar temperature sensitivity to the passive case.The SpongeBOB model is significantly less sensitive to temperature than the Davie and Buffett (2001) onedimensional model (Fig. 17a).These are essentially equilibrium calculations, comparable to those of Buffett and Archer (2004).
The difference is probably due to the fundamentally different processes controlling vertical methane fluxes between the two models.The 1-D model relies on the transport of dissolved-phase methane into the stability zone.In a warmer ocean, the stability zone gets thinner, which increases the diffusive loss of dissolved methane in that model by requiring the diffusive gradient to be steeper.This had the impact in the Buffett and Archer (2004) simulations of restricting the formation of hydrate to regions of higher POC content, reducing the ocean's ultimate equilibrium inventory of methane nearly completely with 3 • C of warming.Methane trapping in the stability zone, when it is transported by bubbles in SpongeBOB, would also be affected by a change in the stability zone thickness, but the mechanism is different and so it is not surprising that the two models would have different temperature sensitivity.
The model also predicts a very strong hydrate dependence on the value of the geothermal heat flux (Fig. 17b), with more hydrate in lower heat flux regimes.A sensitivity as strong as this ought to be detectable in observational data.

Sensitivity to POC
The sensitivity of the hydrate inventory to POC can be seen in the model response to changes in the labile fraction of the sedimenting POC (Fig. 17c).In practice, the respiration zone is thin enough that the POC is not often limited by depletion of the labile POC fraction, so the labile fraction acts like a scaling of an apparent bulk POC degradation rate constant.With the tight correlation between organic carbon age and reaction rate (Middelburg, 1989), this range in effective POC degradation rate constants might optimistically be viewed as a range in the uncertainty in the POC degradation rates.
The model sensitivity to POC deposition was also gauged using time-invariant forcing, by varying the "ocean oxic state" variable through the spectrum from " Atlantic" (low POC concentrations in surface sediments, state = 1) to " Pacific" (2), " OMZ" (3), and " Anoxic" (4).The range in sediment surface POC concentrations is intended to bracket something of the range in the ocean, and Figure 17d shows the sensitivity to that range.Changes in POC deposition have a much greater capacity to increase the hydrate inventory than changes in temperature could, because the ocean is near freezing and it can mostly warm from here.
The ability of high POC to overcome warm temperatures to produce hydrates in the geologic past (Gu et al., 2011) is not supported by the time-dependent and time-invariant simulations shown in Fig. 16c.The simulation without changes in temperature or POC begins to accumulate hydrate with 50 million years of the simulation (dashed line, Fig. 16c).The time-varying simulation, subjected to warmer temperatures but higher POC rain during the middle of the simulation (Fig. 4), does not produce hydrate until the warm, high-POC time is finished at around 150 million years (Fig. 16c, solid line).
The hydrate inventory in the model responds to methane sources other than biological degradation of POC.A model sensitivity run with thermogenic methane production disabled (called No Thermogen), and the simulation No Petroleum, which disables secondary respiration of migrating petroleum, each produce about 60 % of the hydrate from the comparable low-resolution base case.
However, the hydrate inventory is so responsive to reasonable changes in primary POC respiration that were we to neglect these processes in the model, the rate of primary respiration could be increased enough within the uncertainties to bring the hydrate inventory up to much higher values (e.g. the Labile 100 % simulation).Therefore, these results do not constitute a hard prediction of the relative contributions of these processes to the hydrates in real margin sediments, but rather a statement of the relative sensitivities of the hydrate inventory, which is much higher to changes in biological respiration of POC than to petroleum or thermogenic methane.

Sensitivity to bubble transport
The hydrate inventory in the model depends critically on methane transport by rising bubbles through the sediment column, and re-dissolution if the bubbles enter an undersaturated condition.This mechanism is novel to SpongeBOB, not shared by the Davie and Buffett (2001) or most other 1-D methane hydrate accumulation models.If bubble transport is disabled, hydrate nearly disappears from the results (Fig. 16, simulation No Bubb Mig). Figure 17e shows the extreme model sensitivity to the re-dissolution height scale, which determines what fraction of the methane will be recaptured by re-dissolution rather than escape directly to the sea floor.When a scale height of 2 km is imposed, hydrate nearly disappears from the simulation.Shortening the height scale from 500 m to 100 m increases the hydrate inventory by a more than a factor of four.Clearly this is a parameter to which the model is very sensitive, and further work will be required to capture the complexities of the real world in a numerical model.

Pore fluid flow heterogeneity
Old iodine-129 ages in pore waters at Blake Ridge, relative to the age of the solids, can only be reproduced in the model as resulting from heterogeneity in the pore fluid flow field, as crudely mimicked in the code by the vertical high-permeability channels.This conclusion is reinforced by the heterogeneity in penetration depth for dissolved sulfate, SO 2− 4 (Borowski et al., 1996), which seem to show a bimodal distribution, with some sites where the methane/sulfate boundary is just a few tens of meters, and others where sulfate reaches hundreds of meters or more.The shallow sites correlate with the presence of methane hydrate, just like in the model where the hydrate is found.The importance of the channels justifies the imposition of an upward flow in the Davie and Buffett one-dimensional model ( 2001), which was interpreted as representing a zone of upward flow, analogous to the channels here.
However, the impact of the channels on the hydrate inventory itself is relatively small in the model and not apparently robust, in that the high resolution model saw an increase in hydrates when the channels were disabled, but the lower resolution simulations showed no difference.Disabling the horizontal anisotrophy in sediment permeability along with the permeable channels (No Anisotrophy) also had a minor impact.

Other model sensitivities
The Low Slope scenario produces significantly more hydrate than the base case, probably because of the expanded seafloor area within the hydrate-bearing zone on the broader slope of that simulation.Eliminating landslides (No Slides), which allowed a somewhat steeper slope but did not change the hypsometry very much, did not change the hydrate inventory relative to the base case.
The No SO4 scenario increased the hydrate inventory by about 70 %, by steering respiration toward methanogenesis rather than sulfate reduction, and eliminating anaerobic oxidation of methane as a methane sink availability of dissolved sulfate in the ocean as an upper boundary condition to the model.
Most of the other scenarios had a fairly minor impact on the hydrate inventory, when viewed within the context of the limitations of the coarse numerical grid.Examples include clay dewatering and dissolved organic carbon production and degradation.

Conclusions
The SpongeBOB model treats the coastal margin methane cycle from a much larger perspective than previous models of methane hydrate accumulation in margin sediments.The model domain spans to bedrock, closing the budgets for pore fluid and carbon, and into the onshore-offshore dimension.When heterogeneities in the sediment permeability are introduced, the lateral dimension becomes an integral part of the dynamics, rather than the model being a sort of palisade of one-dimensional systems.Also, rates of sedimentation are governed by the evolution of the water depth overlying the sediment column, so that sedimentation variations in time are controlled by things going on in the lateral dimension.Its two-dimensional formulation therefore allows the model to span geologic time.
Bubble migration through the sediment column is an important part of the methane cycle in the model.Viewed from a basin perspective, it is no surprise that methane is lost from the deep sediment column.Another related but separate question is how much of the methane gas lost from the deep sediment column is recaptured to the dissolved phase.
The SpongeBOB model departs from most previous models of methane hydrate accumulation (Burdige, 2011;Chatterjee et al., 2011;Davie and Buffett, 2001;Davie and Buffett, 2003a;Davie and Buffett, 2003b;Garg et al., 2008;He et al., 2011;Liu and Flemings, 2007;Xu, 2004) by postulating that this upward bubble flux might play a role in the stratigraphic methane hydrate deposits on passive continental margins.It turns out that the model methane hydrate inventory depends very strongly on the scale height for methane re-dissolution within an undersaturated zone.The fate of methane gas within the sediment column clearly needs further attention and development.
Heterogeneity in the vertical flow field, simulated very crudely in SpongeBOB as vertical channels of increased permeability, affect chemical tracers such as iodine-129 and SO 2− 4 , as well as focusing the formation of methane hydrate in the channels.Old I 129 ages in hydrate-bearing sediments (Fehn et al., 2000) and the observed correlation between SO 2− 4 penetration depth and presence of hydrate (Borowski et al., 1996) are consistent with the conclusion that hydrate distribution is steered by upward flow.
The abundance of methane hydrate in the model is poised on a knife's edge, responding very strongly to changes in the degradation rate of organic matter, to gas transport and redissolution in the sediment column, and to ocean and subsurface temperatures.For this reason it is currently impossible to calculate the global inventory of methane hydrate from modeling first principles.The goal of the process-oriented model

Fig. 2 .
Fig. 2. Sediment transport model results.(A) Mean grain size (φ = log 10 (mm)) in the sediment column half-way through (left) and at the end of the simulation (right).(B) Snapshots of the sedimentation rate of particles of various sizes.(C) Snapshots of the total sedimentation rate, with a depocenter just off the shelf break.An animation of this figure is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/fig2.atlantic.movie.gif.
Fig. 3. (A) Sediment column constructed from a seismic section off Cape Cod, from Kennett (1982), with cutout (B) resized for comparison with model results (C).

Fig. 4 .
Fig. 4. (A)Sea level time variation with time imposed on the simulation, and corollary changes in the "ocean state", which drive POC deposition concentration patterns as in (B), where Atlantic and Pacific represent conditions typical of the East or Northwest coasts of United States, "OMZ" is stronger oxygen minimum conditions typical of the Eastern equatorial Pacific, and Anoxic refers to anoxic conditions.Concentrations of POC in ocean sediments vary through this range and beyond; the intent here is to allow for variations that span something like the observed range.

Fig. 5 .
Fig. 5. Porosity at the end of the simulation, section (A) and profile (B) taken on the heavy solid line in (A): (solid line) full model solution, (long dashes) "drained" porosity contribution from β 1 , (short grey dashes) sum of drained porosity contribution from β 1 + β 2 .

Fig. 6 .
Fig.6.Log of the vertical permeability from the end of the simulation, (A) from the Base case and (B) the No Chan case with the vertical permeable channels disabled and at lower horizontal resolution.

Fig. 7 .
Fig. 7. Vertical velocity relative to the moving sea floor (w seafloor ) for the Base (A) and No Channels (B) cases, from the end of the simulation.An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/w seafloor.movie.gif.

Fig. 8 .
Fig. 8.Rates of vertical flow at the sea floor.Solid lines are the Base scenario, dashed the No Channel run.(A) Darcy flows (w Darcy ) are upward, while (B) flows relative to the moving seafloor (w seafloor )can be found in both directions.An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/w seafloor sf.movie.gif.

Fig. 9 .
Fig.9.POC concentrations (dry weight percent) (A and C) and H/C elemental ratio of organic matter (B and D).From half-way through (A and B) the simulation, at a period of high sea level, and (C and D) at the end of the simulation, when sea level is relatively low, as today.An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/poc.movie.gif.
et al.: A model of passive margin carbon and methane function of the true Heaviside function.

Fig. 10 .
Fig. 10.Methane sources in the sediment column.(A) From biological respiration of POC.(B) From biological degradation of migrating petroleum.(C) From thermal alteration of POC.An animation of the simulation can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/ch4 src.movie.gif.

Fig. 11 .
Fig. 11.Concentration of dissolved methane plotted in moles m −3 (A and B), and relative to saturation (C and D).Base model run (A and C), and No Channels (Band D).An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/spongebob passive/ch4.movie.gif.

Fig. 12 .
Fig. 12. Bubbles (A, C, E) and hydrate (B, D, F) (volume percent of pore space).From the Base model (A and B), the No Channels model (C and D), and with bubble migration disabled (E and F).An animation of the simulation is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/bubbles.movie.gif.

Fig. 14 .
Fig.14.Carbon isotopic composition of methane from half-way through (A), and at the end of the simulation (C).(B) is the δ 13 C of DIC at the end of the simulation, compared with data from Blake Ridge ODP sites 994 and 997(Rodgiruez et al., 2000) and the West African margin Sivan et al. (2007) (D).The model lacks the resolution to capture the isotopically light spike just below the sea floor, but the deep isotopic composition matches well.Animations of this plot is archived in the Supplement and can be viewed at http: //geosci.uchicago.edu/∼ archer/spongebob passive/del13.movie.gif.
Fig. 15.Iodine-129 porewater age -sediment age, in millions of years.With permeable channels 100 100 Myr of model time, (A) and 200 Myr (B), without permeable channels (C), and (D) with pore fluid flow disabled, t = 30 myr instead of 200 myr for the other simulations.Animations of this plot is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/i129.movie.gif.

Fig. 16 .
Fig. 16.Results of model sensitivity runs.(A and B) Hydrate inventory at the end of the simulation, (A) for the time-uniform forcing runs and (B) for the time-varying forcing runs.(C) The time trajectory of hydrate inventory, solid line is Base (time-varying), dashed line is simulation T 0 (time-uniform).

Fig. 17 .
Fig. 17.Summary of model sensitivities.Y-axis is the hydrate inventory relative to the appropriate base case = 1.(A) Temperature, (B) geothermal heat flux, (C) percent labile fraction of depositing POC, (D)ocean oxic state (see Fig. 3), and (E) migrating bubble re-dissolution scale height.

Table 1 .
Model scenarios with uniform sea level and POC deposition state (ocean oxic state = 3).

Table 2 .
Simulations with a time-varying sea level and ocean oxic state (POC deposition).

Archer et al.: A model of passive margin carbon and methane cell
to another as the vertical grid in the model stretches to fill the expanding sediment column (explained above).The fluid component of this bookkeeping "flow" is calculated as the change in elevation of a grid boundary in a time step, multiplied by the porosity.
-129 porewater age -sediment age, in millions of years.With permeable channels 100 100 Myr of model time, (A) and 200 Myr (B), without permeable channels (C), and (D) with pore fluid flow disabled, t = 30 myr instead of 200 myr for the other simulations.Animations of this plot is archived in the Supplement and can be viewed at http://geosci.uchicago.edu/∼ archer/ spongebob passive/i129.movie.gif.