Methane and sulfate dynamics in sediments from mangrove-dominated tropical coastal lagoons, Yucatán, Mexico

Abstract

Porewater profiles in sediment cores from mangrove-dominated coastal lagoons (Celestún and Chelem) on the Yucatán Peninsula, Mexico, reveal the widespread coexistence of dissolved methane and sulfate. This observation is interesting since dissolved methane in porewaters is typically oxidized anaerobically by sulfate. To explain the observations we used a numerical transport-reaction model that was constrained by the field observations. The model suggests that methane in the upper sediments is produced in the sulfate reduction zone at rates ranging between 0.012 and 31 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$, concurrent with sulfate reduction rates between 1.1 and 24 mmol SO${}_{\mathrm{4}}^{\mathrm{2}-}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$. These processes are supported by high organic matter content in the sediment and the use of non-competitive substrates by methanogenic microorganisms. Indeed sediment slurry incubation experiments show that non-competitive substrates such as trimethylamine (TMA) and methanol can be utilized for microbial methanogenesis at the study sites. The model also indicates that a significant fraction of methane is transported to the sulfate reduction zone from deeper zones within the sedimentary column by rising bubbles and gas dissolution. The shallow depths of methane production and the fast rising methane gas bubbles reduce the likelihood for oxidation, thereby allowing a large fraction of the methane formed in the sediments to escape to the overlying water column.

Introduction

Wetlands are the largest natural source of methane (CH${}_{\mathrm{4}})$ to the atmosphere, accounting for between 20 and 25 % of the global atmospheric methane budget (Fung et al., 1991; Whalen, 2005). Methane produced in wetlands is primarily biogenic, originating from microbial activity in anaerobic sediments and soil. Since sulfate-reducing bacteria outcompete methanogens for common substrates (Oremland and Polcin, 1982), freshwater wetlands typically have much higher methane fluxes to the atmosphere than brackish to fully marine wetlands (Bartlett et al., 1985, 1987; Segarra et al., 2013). Marine and estuarine sediments are generally characterized by comparatively low rates of methanogenesis with a methane production and accumulation zone located deeper within the sediment below the sulfate reduction zone (Holmer and Kristensen, 1994; Martens and Klump, 1984; Poulton et al., 2005; Segarra et al., 2013). In these marine or estuarine systems methane that diffuses upwards towards the sediment surface can be oxidized both anaerobically (AOM) and aerobically within the sediments and in the water column, reducing emissions to the atmosphere (Whalen, 2005).

Maps of (a) the Yucatán Peninsula with lagoon locations, (b) Celestún Lagoon and (c) Chelem Lagoon showing the sampling stations (circles) of sediment cores.

[Figure omitted. See PDF]

Despite brackish to marine salinities, methane fluxes comparable to those measured in freshwater wetlands have been reported for coastal mangrove-dominated lagoon systems in several places around the world, including Florida (Barber et al., 1988), Puerto Rico (Sotomayor et al., 1994), India (Biswas et al., 2004, 2007; Purvaja and Ramesh, 2000, 2001; Ramesh et al., 1997, 2007; Verma et al., 1999), Tanzania (Kristensen et al., 2008), Thailand (Lekphet et al., 2005), China (Alongi et al., 2005), the Andaman Islands (Linto et al., 2014) and Australia (Call et al., 2015). The anaerobic and organic-rich sediments found in these systems provide a suitable environment for methanogenesis, yet the extensive supply of sulfate from seawater should favor sulfate reducers over methanogens in the shallow sections of the sediments (Kristensen et al., 2008; Lee et al., 2008). There are, however, several possible ways for coastal mangrove lagoons to sustain relatively high methane fluxes despite high sulfate concentrations. For example, if the microbial activity of sulfate reducers is high and sulfate replenishment from the overlying water is slow, sulfate may become depleted in the upper centimeters of the sediment, thus allowing methanogenesis to occur close to the sediment surface. Additionally, methanogens can co-exist with sulfate reducers when non-competitive substrates (those used only by methanogens and not by sulfate reducers) are available. Moreover, in some systems methane may migrate from deeper in the sediment to shallower depth and to the water column. Typically, a large percentage of the methane produced in sediments is oxidized prior to reaching the atmosphere, and in shallow-water systems, the oxidation takes place primarily in the sediments and not in the water column (Martens and Klump, 1980; Mitsch and Gosselink, 2000; Weston et al., 2011; Segarra et al., 2013, 2015). However, accumulation and transport of methane in gas bubbles reduces the exposure time of methane to oxidants such as oxygen and sulfate, allowing a large fraction of gas to escape the sediment (Barnes et al., 2006; Martens and Klump, 1980).

The objective of this study was to examine porewater methane distributions within the sediments of two mangrove-dominated coastal lagoons in Mexico and relate them to sulfate concentrations in sediment porewaters. We aim to gain a better understanding of the factors controlling the methane flux from coastal mangrove-dominated lagoon sediments. To this end, we applied a numerical transport-reaction model based on Wallmann et al. (2006) and Chuang et al. (2013) to simulate porewater methane and sulfate concentration profiles. We also performed sediment slurry incubation experiments to test the effect of competitive and non-competitive substrates on methanogenesis in the lagoon sediments. The results provide quantitative data on methane dynamics in coastal mangrove-dominated lagoon systems and highlight their importance as methane sources to the atmosphere.

Study sites

Fieldwork was conducted in two mangrove-dominated coastal lagoons located on the western Yucatán Peninsula, Mexico (Fig. 1). The typical climatological pattern for this area consists of a dry season (March–May), a rainy season (June–October) during which the majority of the annual rainfall (> 500 mm) occurs, and the “nortes” season (November–February), which is characterized by moderate rainfall (20–60 mm) and intermittent high wind speeds greater than 80 km h${}^{-\mathrm{1}}$ (Herrera-Silveira, 1994).

Celestún Lagoon (20${}^{\circ }$52${}^{\prime }$ N, 90${}^{\circ }$22${}^{\prime }$ W) is long, narrow, and relatively shallow (average depth $=$ 1.2 m). The inner and middle sections of the lagoon always have lower salinities than the section near the mouth due to year-round discharge of brackish groundwater from multiple submarine springs (Young et al., 2008). Salinity within the lagoon fluctuates seasonally, with salinity in the inner zone ranging from 8.9 to 18.2 during the course of this study, grading out to marine salinities at the mouth of the lagoon (Young et al., 2008). The lagoon is surrounded by 22.3 km${}^{\mathrm{2}}$ of a well-developed mangrove forest, and has experienced relatively little disturbance from human development and/or pollution such as wastewater discharge (Herrera-Silveira et al., 1998). Sediments in Celestún consist primarily of autochthonous carbonate ooze.

Chelem Lagoon (21${}^{\circ }$15${}^{\prime }$ N, 89${}^{\circ }$45 W; average depth $=$ 0.7 m), in contrast, receives very little groundwater input and the surrounding area has been heavily impacted by urban development. Salinity in Chelem ranged from brackish to hypersaline (24.8–40.3 during the study period), and vegetation surrounding the lagoon consists of scrub mangrove forest (Herrera-Silveira et al., 1998). The construction of Yucalpeten Harbor in 1969 (Valdes and Real, 1998) increased the circulation and resulted in sandy marine sediments entering the lagoon. Sediments in Chelem deposited since 1969 consist of a heavily bioturbated mix of sands and autochthonous carbonate ooze, with a large number of shells of living and dead burrowing organisms (Valdes and Real, 1998). In the following text, CEL and CH denote Celestún Lagoon and Chelem Lagoon, respectively.

Sampling and analytical methods

Porewater solutes

Sediment cores were collected along lengthwise transects in both lagoons during the three different seasons; April 2000 (dry season), December 2000 (nortes season), and October 2001 (late rainy season). Duplicate samples (1_1CH_Oct01 and 1_2CH_Oct01) were collected at station 1CH in Chelem lagoon. Sediments were sampled using hand-held acrylic push cores (7 cm inner diameter) either 30 or 60 cm in length. The push cores had holes drilled along the side at 2 cm intervals, which were sealed with electrical tape prior to sampling. Subsamples for porewater methane analysis were collected in the field immediately after core collection from the holes along the sides of the push cores, using plastic 3 mL syringes with the needle attachment end removed. The sediment plugs from the syringes were immediately extruded into 20 mL glass Wheaton bottles and sealed with blue butyl stoppers and aluminum crimp caps. 3 mL of degassed Milli-Q water and 0.3 mL of saturated mercuric chloride (HgCl${}_{\mathrm{2}})$ solution were added to create a slurry and halt all biological activity within the sample.

After subsampling, the cores were capped, the holes were resealed, and the cores were transported back to the lab for sectioning and porewater extraction. The cores were extruded and sliced into 2.5 cm depth intervals in an anaerobic glove bag under an N${}_{\mathrm{2}}$ atmosphere and transferred into centrifuge tubes for porewater extraction. Core length was measured immediately after collection and just prior to extrusion in order to correct for compaction during transport. Average compaction was 6 % of the total core length, and never exceeded 20 %.

Porewater for sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-})$ and chloride (Cl${}^{-})$ analyses was extracted by centrifuging all the sediment from each depth interval and filtering the porewater through sterile 0.20 $\mathrm{\mu }$m syringe filters. Samples were kept frozen in 20 mL acid-cleaned glass scintillation vials until analysis. Porewater sulfate and chloride concentrations were measured by ion chromatography using a Dionex DX-500 IC equipped with an Ionpac AS9-HC column (4 mm) and AG9-HC (4 mm) guard column. The samples were diluted 5-fold with Milli-Q water prior to analysis in order to bring the sulfate and chloride within the appropriate analytical range for the ion chromatograph.

Methane concentrations for all samples were measured on an SRI 310 Gas Chromatograph (GC) equipped with a flame ionization detector and an Alltech Haysep S 100/120 column (6${}^{\prime }$ $\times$ 1/8${}^{\prime \prime }$ $\times$ 0.085${}^{\prime \prime }$). Helium was used as the carrier gas at a flow rate of 15 mL min${}^{-\mathrm{1}}$ and the column and detector temperatures were maintained at 50 and 150 ${}^{\circ }$C, respectively. Peak integration was performed using Peak Simple NT software. Methane gas standards were prepared by diluting 100 % methane in helium, and five standards bracketing the range of sample concentrations were measured at the beginning, middle, and end of each set of analyses. Average standard error of repeat injections of standards throughout a sample run (between 2 to 6 h of continuous analysis) was 1.8 % ($n=$ 152). Porewater methane concentration in the sediment core subsamples was determined after vigorously shaking the sealed serum bottles containing the sediment slurries to ensure complete mixing, followed by at least 3 min of standing equilibration time to ensure that the porewater methane was fully equilibrated with the headspace in the serum bottles. A small volume of headspace (0.25–0.5 $\mathrm{\mu }$L) was drawn out of each serum bottle using a gas-tight syringe, and analyzed for methane concentration on the SRI 310 GC. The total volume of porewater in each sample was calculated using the difference between the total wet weight of the sediment minus the dry weight of the sediment, correcting for the added water and HgCl${}_{\mathrm{2}}$ solution.

Experimental conditions and sampling time intervals for methane headspace concentration analyses of sediment slurry incubations.

	Treatment	Initial concentration	Experiment	Number of	Methane production
		of treatment	length	measurements	rate (nmol CH${}_{\mathrm{4}}$ cm${}^{-\mathrm{3}}$
			(days)		slurry d${}^{-\mathrm{1}})$
Controls	No amendment (anaerobic)	N${}_{\mathrm{2}}$ headspace	29	3	1.3 $\times$ 10${}^{-\mathrm{4}}$ to 2.0 $\times$ 10${}^{-\mathrm{3}}$
	Autoclaved	N${}_{\mathrm{2}}$ headspace	29	3	0 to 2.6 $\times$ 10${}^{-\mathrm{3}}$
	Aerobic-O${}_{\mathrm{2}}$ gas	16 % O${}_{\mathrm{2}}$ headspace (0.36 mM)	29	3	5.7 $\times$ 10${}^{-\mathrm{4}}$ to 3.5 $\times$ 10${}^{-\mathrm{3}}$
	BES	40 mM	29	3	0 to 1.3 $\times$ 10${}^{-\mathrm{4}}$
Competitive	H${}_{\mathrm{2}}$ gas	100 % headspace (1.8 mM)	29	3	5.4 $\times$ 10${}^{-\mathrm{3}}$ to 6.2
substrates	Acetate	10 mM	29	3	6.8 $\times$ 10${}^{-\mathrm{4}}$ to 9.2 $\times$ 10${}^{-\mathrm{2}}$
	Formate	10 mM	29	3	6.9 $\times$ 10${}^{-\mathrm{4}}$ to 1.6 $\times$ 10${}^{-\mathrm{1}}$
Noncompetitive	Methanol	10 mM	29	4	2.0 $\times$ 10${}^{-\mathrm{2}}$ to 19
substrates	TMA	10 mM	29	4	5.4 $\times$ 10${}^{-\mathrm{4}}$ to 40

Depth profiles of modeled (lines), measured (circles) and calculated (triangles) concentration of dissolved methane (dashed lines; open circles), sulfate (solid lines; solid circles) in the upper panel and sulfate depletion (solid lines; solid triangles), zero sulfate depletion (dashed lines) and chloride concentration (gray circles) in the lower panel for each profile type (Groups 1–4, see text). One selected profile per group is shown here for illustration and the other profiles for each group (9 cores for Group-1, 6 cores for Group-2, 2 cores for Group-3 and 3 cores for Group-4) are presented in the Appendix A (Fig. A1). CEL and CH represent cores collected from Celestún Lagoon and Chelem Lagoon.

[Figure omitted. See PDF]

Sediment slurry incubation experiments

Sediment slurry incubations were performed in order to examine changes in methane production over different time intervals and at different substrate concentrations (Table 1). Incubations consisted of three competitive substrates (H${}_{\mathrm{2}}$, acetate, formate), two non-competitive substrates (methanol, trimethylamine (TMA)), and four types of controls. The controls (preparation methods are described below) consisted of an un-amended sediment control under anaerobic conditions, an un-amended aerobic control (partial oxygen headspace), a killed control in which the sediment was autoclaved to kill all living organisms in the sediment, and a chemical control in which biological methanogenesis was inhibited through the addition of 2-bromoethanesulfonic acid (BES) to a final concentration of 40 mM within the slurry. Triplicate bottles were prepared for each condition (controls and substrate additions), and methane headspace concentrations were measured at 3–4 time intervals over the course of 29 days.

All the sediment slurries were prepared semi-anaerobically by homogenizing the sediment in a blender with an artificial seawater mixture in a $\mathrm{1}:\mathrm{1}$ ratio under continuous flow of nitrogen gas. Large pieces of leaves, twigs, and shells were removed from the sediment prior to homogenization. 70 mL glass Wheaton bottles were flushed with nitrogen gas for 1 min prior to the addition of the sediment slurry. 30 mL of slurry was then added to each bottle under continuous nitrogen flow, and the bottles were sealed using blue butyl rubber stoppers and aluminum crimp seals. Substrate additions were made by injecting the substrate solution into the bottle immediately after sealing the bottles, except for the H${}_{\mathrm{2}}$ gas treatment and the aerobic control. For the addition of H${}_{\mathrm{2}}$, the entire headspace of the bottles was flushed with 100 % H${}_{\mathrm{2}}$ gas. After each headspace sampling the H${}_{\mathrm{2}}$ gas removed by microbial activity in the sediment was replaced by inserting a gas tight syringe filled with 100 % H${}_{\mathrm{2}}$ gas into the bottles, and allowing the gas to be drawn into the bottles until equilibrium pressure was reached. The aerobic controls were prepared like the anaerobic un-amended controls, except that 8 mL (20 % of the total headspace) of 100 % O${}_{\mathrm{2}}$ was added to the bottles immediately after they were sealed. In order to ensure that the sediment slurries remained aerobic, 100 % O${}_{\mathrm{2}}$ was added to the bottles throughout the incubation period. The sediment slurries were kept at room temperate (22 ${}^{\circ }$C) and agitated continuously on a shaker table throughout the course of the incubations.

Headspace samples (0.25 mL) were extracted from the bottles at each time interval using a gas-tight syringe. Methane concentrations were measured on an HP 5730A GC equipped with a flame ionization detector. GC calibration and creation of standard curves were based on successive dilutions of 100 % methane. Analytical error was approximately 5 % for methane concentrations below 10 ppm-v (446 nM), and less than 3 % for methane concentrations above 10 ppm-v as determined by repeat analyses of standards and samples.

Results

Porewater concentrations of dissolved species

Representative porewater methane profiles were plotted alongside sulfate profiles in Figs. 2 and A1 in Appendix A. Profiles were assigned to one of four profile types based on the relation between methane and sulfate distributions down core (see below). Considerable spatial and temporal variability in porewater chemistry was observed with no systematic seasonal differences in concentration trends. For example, porewater methane concentrations varied by up to 3 orders of magnitude in both lagoons, even between sites in close proximity to each other (i.e. 1CEL and 2CEL, Oct01; 1CH and 2CH, Dec00), and at the same station sampled during different seasons (i.e. 2CEL Dec00, Oct01; 1CH Apr00 and Oct01). No consistent differences were evident between the stations at the sides of the lagoons and those located in the center of the lagoons, or between stations located in the inner zone of the lagoons and those located near the mouth. For instance, supersaturated methane concentrations of 1.1 and 1.3 mM were observed in cores 1CEL_Jul02 (the inner zone of Celestún lagoon) and 14CEL_Dec00 (near the mouth of Celestún lagoon), respectively. This is particularly interesting because the water column at the mouth of the lagoon has much higher salinities than the water of the inner zone (Young et al., 2008). The variability (both spatial and temporal) in the porewater methane concentrations and in profile types suggests a very dynamic system where both concentration and distribution patterns in the porewater vary constantly. Such variability is indicative of rapid methane production and efflux rates.

Relationship between (a) [Cl${}^{-}$] and [SO${}_{\mathrm{4}}^{\mathrm{2}-}$], (b) [Cl${}^{-}$] and [CH${}_{\mathrm{4}}$] and (c) [CH${}_{\mathrm{4}}$] and [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] in porewater samples.

[Figure omitted. See PDF]

Porewater sulfate concentrations ranged from 0.21 to 35.3 mM in Celestún lagoon and from 4.13 to 33.5 mM in Chelem lagoon and showed different trends (Figs. 2, A1). In many of the cores a negative relation between methane and sulfate was observed. Specifically, higher sulfate was associated with lower methane in cores located near the mouth of the lagoons (16CEL_Jul02, 16CEL_Oct01, 14CEL_Oct01, 14CEL_Jul02 and 5CH_Apr00) and lower sulfate with high methane in the inner zone of the lagoons (e.g. cores 1CEL_Jul02, 1CEL_Dec00, 3CEL_Jul02, 3CEL_Apr00, 1_1CH_Oct01, and 1_2 CH_Oct01).

The relationship between porewater salinity (represented by chloride concentration), methane and sulfate concentrations was spatially and temporally variable (Fig. 3). Generally, higher sulfate concentrations were associated with higher chloride in cores located near the mouth of the lagoons and lower sulfate with lower chloride in the inner zone of the lagoons (Fig. 3a). Despite these general trends there were no clear consistent relationships between methane and chloride (Fig. 3b) and sulfate and methane (Fig. 3c) when the data were considered collectively. The lack of consistent trends suggests multiple processes impacting the distribution of methane and sulfate. These include physical processes, such as mixing and dilution by seawater or groundwater, and biological processes such as sulfate reduction, methanogenesis and methane oxidation. Brackish groundwater enters Celestún lagoon through at least 30 subsurface discharge points (Young et al., 2008), and the chloride profiles suggest that some of this groundwater may seep through the sediments, resulting in localized decline in porewater salinities.

To account for mixing with seawater or freshwater and to extract information on the biological and chemical processes controlling the distribution of porewater solutes, the observed sulfate depletion ([SO${}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}$]${}_{\text{OBS}})$ relative to seawater was calculated as the difference between the expected sulfate concentration contributed by seawater (based on porewater chloride concentration) and the measured sulfate concentration: $$[{\mathrm{SO}}_{{\mathrm{4}}_{\mathrm{dep}}}^{\mathrm{2}-}{]}_{\mathrm{OBS}}={\displaystyle \frac{[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}{]}_{\left(\mathrm{SW}\right)}}{[{\mathrm{Cl}}^{-}{]}_{\left(\mathrm{SW}\right)}}}\times [{\mathrm{Cl}}^{-}{]}_{\left(\mathrm{measured}\right)}-[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}{]}_{\left(\mathrm{measured}\right)}$$ where 0.05171 is taken as the $\frac{[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}{]}_{\left(\mathrm{SW}\right)}}{[{\mathrm{Cl}}^{-}{]}_{\left(\mathrm{SW}\right)}}$ ratio (Pilson, 1998). Positive values indicate that sulfate has been removed from the porewater, most likely through sulfate reduction, while negative values indicate an external source of sulfate not associated with chloride, in this case groundwater (see discussion below).

Based on the observed trends in sulfate depletion, when considered together with methane, four different porewater trends can be described, referred to as Groups 1 through 4 (Figs. 2, A1). The majority of profiles fell into Group-1 (ten cores); these profiles showed positive sulfate depletion profiles (e.g. sulfate consumption or loss) with methane profiles mirroring the sulfate concentration profiles (methane production or input). The peaks for methane and sulfate depletion occurred at the same depth as the lowest measured sulfate concentrations. In Group-2 (seven cores), sulfate depletion also showed positive values (sulfate consumption) but not throughout the core. In some cores sulfate depletion was close to zero at shallow depths and then increased with depth and in other cores positive sulfate depletion values appeared at the surface of the sediments and then decreased to almost zero at deeper depths. Methane concentrations for this group showed no clear relation to the sulfate profiles. In Group-3 (three cores), sulfate depletion showed negative values (e.g. sulfate addition). The values became more negative toward the deeper sediment starting from zero right at the surface suggesting that sulfate was being added from the bottom of the sediment section. In Group-4 (four cores), there was almost no sulfate depletion (sulfate concentrations similar to seawater) from the surface to the deeper depths and methane concentrations were low (< 0.25 mM) increasing at depth, indicating a deeper source of methane.

(a) Headspace methane concentrations in sediment slurry incubations. (b) Expansion of (a), showing results for acetate, formate, and controls. (c) Expansion of (a), showing results for controls only. Error bars represent 1 standard deviation for triplicate sample bottles.

[Figure omitted. See PDF]

Sediment slurry incubation experiments

All the sediment slurries with added substrates showed an increase in headspace methane concentration that was significantly greater than those observed with either the un-amended aerobic and anaerobic controls or the treated controls (Fig. 4). The greatest increases in headspace methane concentration were seen with additions of the two non-competitive substrates, TMA and methanol. The H${}_{\mathrm{2}}$ treatment showed the next highest methane production rate, followed by formate and acetate. Of the four control conditions, the un-amended, anaerobic treatment had the highest overall increase in headspace methane concentration. The aerobic treatment had an initial higher increase in headspace methane concentration than the un-amended, anaerobic treatment, although there was no detectable change in the headspace methane concentration in the aerobic treatment between 150 and 700 h. Both the autoclaved and BES treatments did not show any changes in headspace methane concentration greater than the instrumental detection limits. The maximum methane production rates for each treatment are listed in Table 1.

Discussion

Co-existence of methane and sulfate in sediments

Seawater transport into the sediment by diffusion and bioirrigation due to the activity of burrowing animals has clear effects on porewater solutes. These processes are a source of seawater sulfate and mask sulfate loss by microbial reduction. Although, as indicated above, considerable variability in porewater profile distribution trends was observed, and different profile types were found throughout the lagoons, certain trends were more common at distinct locations. Specifically, sites characterized by sulfate addition from input of seawater into the sediment (cores in Group-4) were found primarily near the mouth of both lagoons where low methane was associated with near-zero sulfate depletion. Negative sulfate depletion (Group-3), on the other hand, which indicates the presence of porewater that is enriched in sulfate relative to chloride, was seen primarily in the middle zone of Celestún Lagoon where groundwater springs rich in sulfate due to anhydrite dissolution are present, as reported by Perry et al. (2002, 2009). Positive sulfate depletion profiles co-occurring with methane (Groups 1 and 2) were seen throughout the lagoons but mostly at sites in the inner zone of both lagoons, suggesting significant sulfate reduction at rates higher than the replenishment from sulfate rich groundwater or from the overlying seawater and a source of methane to the shallow sections of the sediment.

It is surprising that at many sites, particularly within Groups 1 and 2 in the inner zone of both lagoons (1CEL, 2CEL, 3CEL and 1CH), high concentrations of methane and sulfate co-occurred at the same depth in the sediment. Co-existence of methanogenesis and sulfate reduction is not normally observed because sulfate reduction is more energetically favorable than methanogenesis, and sulfate reducers should outcompete methanogens for common substrates such as hydrogen and acetate (Oremland and Polcin, 1982; Jørgensen and Kasten, 2006). Moreover, anaerobic oxidation of methane (AOM) coupled with sulfate reduction at the base of the sulfate reducing zone should further deplete methane (Capone and Kiene, 1988; Valentine and Reeburgh, 2000). There are several possible explanations for these observations: firstly, the high methane concentrations measured in the sulfate rich porewater may be supplied by a rapid non-diffusive mechanism from below the sulfate reduction zone (like rising gas bubbles), limiting the exposure time to AOM. Secondly, methane may be produced in situ at these depths supported by a high abundance of competitive substrates in the sulfate reduction zone hence sustaining both methanogenesis and sulfate reduction (Holmer and Kristensen, 1994). Thirdly, methanogens may instead be able to thrive on various non-competitive substrates (Oremland and Polcin, 1982; Wellsbury and Parkes, 2000; Lee et al., 2008; Taketani et al., 2010). Indeed, use of non-competitive substrates by methanogens, including methanol, trimethylamines and dimethylsulfide, has been reported for mangrove sediments, coastal lagoons and continental shelf sediments (Ferdelman et al., 1997; Lyimo et al., 2000; Mohanraju et al., 1997; Purvaja and Ramesh, 2001; Torres-Alvarado et al., 2013; Maltby et al., 2016). Our slurry incubation experiments demonstrated that the methanogenic community at Celestún is capable of using a wide range of substrates, including H${}_{\mathrm{2}}$, acetate, formate, methanol, and trimethylamine (Fig. 4). Both methanol and trimethylamine are not utilized by sulfate reducers, which could allow methanogens to thrive in the sulfate reduction zone (Fig. 4). The use of non-competitive substrates by the methanogenic community has important implications for methane fluxes to the atmosphere as it allows for methane production at shallow depths in the sediment and reduces the potential for complete oxidation of methane. Although processes and trends similar to those described above have been reported for other mangrove sediments (e.g., Lee et al., 2008; Purvaja and Ramesh, 2001), the co-occurrence of sulfate and methane and related biogeochemical reactions in these reports remain qualitative in nature. In the following section, we use a transport-reaction model to better quantify the processes controlling methane fluxes from the sediments in these mangrove-dominated tropical coastal lagoons.

Model set-up and application to mangrove-dominated coastal lagoon sediments

In order to understand methane production and consumption and how these processes relate to sulfate dynamics in the lagoon sediments, we used two different approaches to simulate methane and sulfate porewater profiles.

In the first approach, a transport-reaction model was applied to profiles of Group-1 where methane and sulfate co-occur with no indication of groundwater sulfate input and where sulfate reduction surpasses sulfate addition from seawater (Figs. 2; A1). Data in Group-1 have positive net sulfate depletion rates indicative of sulfate reduction. The sulfate depletion is seen within the zone where methane concentrations are high. In these cores the net sulfate depletion rates can be used to derive the minimum methanogenesis rates (see model details in the Appendix A). Reactions considered in this first approach include organic matter degradation via heterotrophic sulfate reduction, methane production via methanogenesis and methane addition from gas bubble dissolution (Haeckel et al., 2004; Chuang et al., 2013).

A second approach (detailed in the Appendix A) was used for simulating the profiles for Group-2, Group-3 and Group-4 which show no positive net sulfate depletion rates when integrated over the core length. These sites are affected by groundwater input or by considerable irrigation and input of seawater. Here, the link between sulfate and methane reactions is less clear and hard to quantify directly.

The following equation was solved to quantify the rates of reaction and transport of dissolved methane and sulfate in the upper 20 cm of the sediments in both approaches (Berner, 1980; Boudreau, 1997): $$\mathrm{\Phi }\cdot {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial \left(\mathrm{\Phi },\cdot ,D_{\text{s}},\cdot ,\frac{\partial C}{\partial x}\right)}{\partial x}}-{\displaystyle \frac{\partial \left(\mathrm{\Phi },\cdot ,v,\cdot ,C\right)}{\partial x}}+\mathrm{\Phi }\cdot R_c,$$ where $x$ is sediment depth, $t$ is time, $\mathrm{\Phi }$ is porosity, $D_{\text{s}}$ is the solute-specific diffusion coefficient in the sediment, $C$ is the concentration of methane or sulfate in the porewater, $v$ is the burial velocity of porewater and $R_C$ is the sum of reactions affecting $C$ (Table A1 in the Appendix A). Solutes were simulated in moles L${}^{-\mathrm{1}}$ of porewater (M). Details of all the reaction terms and parameters and how they were derived for each of the two approaches are given in the Appendix A. The model assumes steady-state conditions to constrain methanogenesis rates at each site. Considering the observed variability in porewater distributions non-steady state simulations would be desirable, yet this would require continuous monitoring of porewater sulfate, methane and chloride concentrations to evaluate temporal changes in sulfate depletion at each site. These time series data are unavailable, hence the modeled “instantaneous” rates bear uncertainties that currently cannot be quantified accurately.

Model sensitivity analysis of methane concentrations for cores in Group-1 to the different processes controlling methane concentrations in porewaters. Black dashed lines denote the standard simulation results: CH${}_{\mathrm{4}}$ production rate $=+R_{\text{MB}}+R_{\text{M}}$. $R_{\text{M}}$ is methanogenesis, $R_{\text{MB}}$ is methane bubble dissolution, $R_{\text{AOM}}$ is anaerobic oxidation of methane and $R_{\text{SD}}$ is net sulfate depletion.

[Figure omitted. See PDF]

Model derived sulfate depletion and sulfate and methane concentrations are shown in Figs. 2 and A1. Modeled porewater data for Group-1 (the most common trend) show that methane generated from organic matter degradation within the upper sediments is a more important methane source than methane diffusing from below and gas bubble dissolution, as further seen in the results of 1CEL_Jul02 and the sensitivity analysis from 2CEL_Jul02 (Fig. 5a). In 1CEL_Jul02, for example, gas dissolution of methane transported from deeper sediments is not necessary at all to achieve a good model fit to the data, and in-situ methanogenesis alone can reproduce methane concentrations similar to the measured data even though methane concentrations are oversaturated (> 1.1 mM (in situ solubility); Fig. 2). In contrast, the modeled methane profile for 2CEL_Jul02 (black dashed line) arguably does require the inclusion of methane from gas dissolution ($R_{\text{MB}}$; Fig. 5a). In Fig. 5a, the gray dashed and solid lines represent only gas dissolution in the methane reaction terms (no methanogenesis within the modeled 20 cm column) using different gas dissolution constants ($k_{\text{MB}}$ values are 0.2 and 0.5 yr${}^{-\mathrm{1}}$ respectively). The model results shown as the gray dashed line simulate the methane concentrations below 10 cm depth, whereas those shown by the gray solid line reproduce methane concentrations in the upper 5 cm, but neither reproduces the data throughout the whole core. Comparing results considering methanogenesis and gas dissolution (black solid line) and methanogenesis only (black dashed line), it is clear that both methanogenesis and some gas dissolution are needed for reproducing the methane distribution observed in core 2CEL_Jul02. This illustrates the complexity of controlling processes and the dynamic nature and resulting temporal variability in methane fluxes at this and the other sites in the lagoons.

Model-derived depth-integrated turnover rates (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$, dissolved methane fluxes to the water column (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ and contributions of methanogenesis to net methane production (%) and heterotrophic sulfate reduction to POC degradation (%). CEL and CH represent cores collected from Celestún Lagoon and Chelem Lagoon.

	Length of model
	column (cm)	$R_{\text{SD}}=R_{\text{SR}}$	$R_{\text{M}}$	$R_{\text{POC}}$	$R_{\text{MB}}$	$F_{\text{methane (top)}}$	$F_{\text{methane (bottom)}}$	$R_{\text{M}}/(R_{\text{M}}+R_{\text{MB}})$	$R_{\text{SR}}/R_{\text{POC}}$	$R_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$	$R_{{\text{CH}}_{\mathrm{4}}}$
Group-1
1CEL_Apr00	20	3.7	0.13	7.7	0.41	0.59	0.06	25 %	97 %
1CEL_Dec00	20	2.2	1.5	7.4	0	0.94	$-$0.60	100 %	59 %
1CEL_Oct01	20	6.2	0.12	13	0.32	0.40	$-$0.04	27 %	98 %
1CEL_Jul02	20	3.6	8.0	23	0	6.0	$-$1.98	100 %	31 %
2CEL_Dec00	20	1.1	0.05	2.3	0.76	0.54	$-$0.27	5.8 %	96 %
2CEL_Jul02	20	11	0.08	22	0.05	0.11	$-$0.02	63 %	99 %
3CEL_Apr00	20	1.3	0.29	3.2	0.24	0.68	0.15	55 %	82 %
3CEL_Jul02	20	7.1	0.25	15	2.2	3.0	0.63	10 %	97 %
1_1CH_Oct01	13.75	24	31	110	0	11	$-$19	100 %	44 %
1_2CH_Oct01	20	3.0	26	58	0	20	$-$5.6	100 %	10 %
Group-2
1CH_Dec00	20					0.52	$-$7.2			4.5	7.8
1CH_Apr00	20					$\sim$ 0	$\sim$ 0			$-$3.2	$\sim$ 0
2CH_Dec00	20					$\sim$ 0	$\sim$ 0			6.9	$\sim$ 0
5CH_Apr00	20					0.012	$\sim$ 0			21	0.013
2CEL_Oct01	20					11	$-$0.01			3.9	11
14CEL_Jul02	20					0.27	$\sim$ 0			3.7	0.27
16CEL_Dec00	20					$-$0.047	0.013			$-$1.8	$-$0.060
Group-3
5CEL_Apr00	10					0.014	$-$0.01			$-$69	0.028
14CEL_Dec00	20					3.4	$-$0.13			10	3.6
14CEL_Oct01	20					0.088	$-$0.01			2.9	0.10
Group-4
16CEL_Jul02	20					0.096	0.02			6.1	0.072
16CEL_Oct01	20					$\sim$ 0	$\sim$ 0			0.83	$\sim$ 0
7CH_Oct01	20					0.13	$\sim$ 0			2.6	0.14
8CH_Dec00	20					$\sim$ 0	$\sim$ 0			0.85	0.012

$R_{\text{SD}}$ is net sulfate depletion (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of SO${}_{\mathrm{4}}^{\mathrm{2}-})$. $R_{\text{SR}}$ is heterotrophic sulfate reduction (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of SO${}_{\mathrm{4}}^{\mathrm{2}-}$). $R_{\text{M}}$ is methanogenesis (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of CH${}_{\mathrm{4}}$). $R_{\text{POC}}$ is total POC mineralization (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of C). $R_{\text{MB}}$ is gas dissolution (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of CH${}_{\mathrm{4}}$). $F_{{\text{methane}}_{\text{(top)}}}$ is the methane flux across the sediment surface (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of CH${}_{\mathrm{4}}$). Negative values in $F_{{\text{methane}}_{\text{(top)}}}$ represent methane flux into the sediments from the water column and vice versa. $F_{{\text{methane}}_{\text{(bottom)}}}$ is the methane flux across the 20 cm lower boundary (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of CH${}_{\mathrm{4}})$. Negative values in $F_{{\text{methane}}_{\text{(bottom)}}}$ represent methane flux to deep sediments and vice versa. $R_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$ is net sulfate input (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of SO${}_{\mathrm{4}}^{\mathrm{2}-})$ and $R_{{\text{CH}}_{\mathrm{4}}}$ is net methane production (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of CH${}_{\mathrm{4}})$ for cores in Group-2 to Group-4. See Appendix A for further model details.

Model derived depth-integrated turnover rates and fluxes

Table 2 lists the calculated depth-integrated turnover rates and fluxes for the individual cores. For profiles in Group-1, methane sources include methanogenesis within the upper 20 cm and/or methane transported from deeper sections (> 20 cm) via bubble transport and dissolution. Methane can be supported fully by methanogenesis without gas bubble dissolution within the modeled upper 20 cm in cores 1CEL_Dec00, 1CEL_Jul02, 1_1CH_Oct01 and 1_2CH_Oct01. Gas bubble transport from deeper sediments and its dissolution contributes more methane than methanogenesis in cores 1CEL_Apr00, 1CEL_Oct01, 2CEL_Dec00 and 3CEL_Jul02.

Methane sinks include fluxes into the water column or methane diffusion into deeper sediments (> 20 cm) and methane oxidation. Our model shows that the major sink for methane is efflux to the water column accounting for over 90 % of methane produced within the upper 20 cm (e.g. 1CEL_Apr00, 1CEL_Oct01, 3CEL_ Apr00 and 3CEL_Jul02). Model-derived methane fluxes to the water column are listed in Table 2 ($F_{\text{methane (top)}})$ and range from 0.012–20 mmol CH${}_{\mathrm{4}}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$. These rates are similar to or up to 2 orders of magnitude larger than fluxes reported for other mangrove lagoon systems in Florida (0.02 mmol CH${}_{\mathrm{4}}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$, Barber et al., 1988; Harriss et al., 1988), Australia (0.03–0.52 mmol CH${}_{\mathrm{4}}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$, Kreuzwieser et al., 2003), and India (5.4–20.3 mmol CH${}_{\mathrm{4}}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$, Purvaja and Ramesh, 2001). Since all methane depth profile types were observed throughout the year with no obvious trends in spatial and temporal distribution (seasons and sampling locations), our results support the idea that methane fluxes in coastal mangrove lagoon systems respond very dynamically to environmental stimuli.

Sulfate sinks include heterotrophic sulfate reduction and AOM, although the model suggests that AOM plays a minor role compared to heterotrophic sulfate reduction. Sulfate reduction ranges from 1.1 to 24 mmol SO${}_{\mathrm{4}}^{\mathrm{2}-}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ and is the major sink for both sulfate and organic carbon in most cores. Sulfate reduction accounts for 2.2 to 48 mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ of total anaerobic carbon respiration, which is in the same range of values listed in Kristensen et al. (2008) for other mangrove sediments.

Mangrove forests are known to be highly productive systems with the capacity to release high concentrations of dissolved organic matter (DOM) to surrounding sediments and porewaters (Kristensen et al., 2008). Tree litter and subsurface root growth provide further significant inputs of organic carbon to mangrove sediments which are unique for this type of system. The rate of organic matter mineralization ($R_{\text{POC}}$; Eq. A6 in the Appendix A) derived from sulfate depletion ranges from 3.2 to 110 mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$. Although our modeling approach for determining degradation rates is not without uncertainty, it is more accurate than rates derived from down-core trends in organic matter content because of temporal variability in accumulation rates in this area (Gonneea et al., 2004). Particulate organic matter will also contain a high amount of refractory carbon that is not easy to quantify and separate from the bulk pool. The derived degradation rates likely represent the more labile particulate components and labile DOM that was not considered (or measured) in this study. The high calculated organic carbon oxidation rates derived here are thus not unexpected since mangrove systems in general (e.g. Dittmar and Lara, 2001; Dittmar et al., 2006; Lee, 1995; Odum and Heald, 1975) and the lagoons in Yucatán in particular are dominated by high concentrations of DOM, a large fraction of which is likely to be labile (Young et al., 2005).

Depth-integrated methane production or consumption rates ($R_{{\text{CH}}_{\mathrm{4}}})$ and net sulfate inputs ($R_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}})$ calculated from Eqs. (A9) and (A10) for cores in Group-2, Group-3 and Group-4 are listed in Table 2. The methane and sulfate net production/consumption rates ranged from $-$0.060 to 11 mmol CH${}_{\mathrm{4}}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ and $-$69 to 21 mmol SO${}_{\mathrm{4}}^{\mathrm{2}-}$ m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ (negative values indicate net sulfate or methane consumptions while positive values indicate production or addition from external sources). Although sulfate depletion values for cores in Group-2 are positive (e.g. net sulfate reduction), sulfate concentrations at some depths of the porewater are relatively high, suggesting continuous sulfate input from deeper within the sediments or from seawater. Cores in Group-3 and Group-4 show negative or zero sulfate depletion that likely results from high rates of sulfate addition from groundwater (Group-3) or seawater (Group-4), thus prohibiting accurate calculation of sulfate reduction and methanogenesis rates. Although, in theory, H${}_{\mathrm{2}}$S oxidation is a possible source for the excess sulfate, we believe that sulfate-rich groundwater input is a more likely source due to correlation between excess sulfate and excess Sr which has been previously described for groundwater in this region (Young et al., 2008). Perry et al. (2002) identified dissolution of evaporites within the freshwater lens at some Yucatán sites as a probable source of excess sulfate in groundwater using the sulfate-to-chloride ratio (100 $\times$ $\frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{Cl}}^{-}\right]}$). Ratios higher than seawater (average seawater is 10.3) are expected where gypsum/anhydrite dissolution occurs (Perry et al., 2002). Another indicator is the Sr $/$ Cl ratio, which is invariably higher in the Yucatan groundwater than in seawater and indicates dissolution of celestite (from evaporite) and/or aragonite (Perry et al., 2002). The region east and south of Lake Chichancanab, referred to as the Evaporite Region by Perry et al. (2002), is characterized by distinctive topography and high sulfate groundwater concentrations (Perry et al., 2002). The groundwater from the Lake Chichancanab area flows northward into the Celestún Estuary which can be recognized by the progressive decrease in the ratio $\frac{{\frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{Cl}}^{-}\right]}}_{\mathrm{groundwater}}}{{\frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{Cl}}^{-}\right]}}_{\mathrm{seawater}}}$ in water from southeast to northwest (Perry et al., 2009). Some groundwater samples with sulfate concentrations as high as 32 mM were reported in Young et al. (2008) and the Sr and sulfur trends for Celestún lagoon (Young et al., 2008) are consistent with our interpretation that gypsum/anhydrite dissolution in groundwater is the source of excess sulfate in the porewater of Group-3 in Celestún lagoon. Due to the impact of groundwater, our sulfate reduction and methanogenesis rates estimated using the model are minimum rates and independent rates of groundwater discharge into each core are needed for obtaining more realistic estimates in these sites.

Maximum model-derived rates of methanogenesis and sulfate reduction for cores in Group-1 and maximum model-derived rates of methane production and sulfate consumption for cores in Group-2, Group-3 and Group-4. CEL and CH represent cores collected from Celestún Lagoon and Chelem Lagoon.

	Max-$R_{\text{M}}$	Max-$R_{\text{SR}}$
	(nmol CH${}_{\mathrm{4}}$ cm${}^{-\mathrm{3}}$ d${}^{-\mathrm{1}})$	(nmol SO${}_{\mathrm{4}}^{\mathrm{2}-}$ cm${}^{-\mathrm{3}}$ d${}^{-\mathrm{1}})$
Group-1
1CEL_Apr00	9.0	304
1CEL_Dec00	116	559
1CEL_Oct01	7.1	740
1CEL_Jul02	564	1425
2CEL_Dec00	4.9	587
2CEL_Jul02	7.4	1323
3CEL_Apr00	20	405
3CEL_Jul02	26	1227
1_1CH_Oct01	2199	1802
1_2CH_Oct01	1959	1476
Group-2
1CH_Dec00	2531	407
1CH_Apr00	1.6	2687
2CH_Dec00	1.1	2835
5CH_Apr00	2.1	8378
2CEL_Oct01	504	715
14CEL_Jul02	19	394
16CEL_Dec00	2.7	330
Group-4
5CEL_Apr00	63	5212
14CEL_Dec00	1517	1756
14CEL_Oct01	23	1007
Group-5
16CEL_Jul02	10	186
16CEL_Oct01	0.08	599
7CH_Oct01	4.1	940
8CH_Dec00	0.57	230

In addition to depth-integrated rates, Table 3 also includes maximum methanogenesis/methane production (Max-$R_{\text{M}})$ and sulfate reduction/consumption (Max-$R_{\text{SR}})$ rates solved by Eq. (2) in the model. It is encouraging that the maximum methane production rates estimated from TMA, methanol and H${}_{\mathrm{2}}$ additions to sediments in the slurry incubations (Table 1) are similar to model-derived Max-$R_{\text{M}}$ at station 16CEL (Table 3), which is the site from which sediments were collected for the slurry incubations. The rates in the TMA, methanol and H${}_{\mathrm{2}}$ treatments from the slurry incubations (Table 1) and in some of our stations are higher than methane production rates from previously reported coastal freshwater and brackish wetland sediments that were measured using radiolabeled acetate and bicarbonate in slurries (Segarra et al., 2013).

Modeled Max-$R_{\text{M}}$ in some cores were 1–2 orders of magnitude higher than rates derived from the sediment slurry incubations (e.g., cores 1CEL_Jul02, 1_1CH_Oct01, 2CEL_Oct01 and 14CEL_Dec00). Although heterotrophic sulfate reduction generally dominates organic matter degradation, Max-$R_{\text{M}}$ values are even higher than the maximum sulfate reduction rates in some cores (1_1CH_Oct01, 1_2CH_Oct01 and 1CH_Dec00). Both the methanogenesis rates measured in the sediment slurry incubations and the modeled maximum methanogenesis rates in this study area were much higher than those reported for some mangrove systems (e.g., Thailand, Kristensen et al., 2000; Malaysia, Alongi et al., 2004; Australia, Kristensen and Alongi, 2006) but similar to other sites in India (Ramesh et al., 2007).

AOM is expected to play an important role in oxidizing methane in tropical porewaters with abundant methane and sulfate (Biswas et al., 2007). However, our model results and sensitivity analyses indicate that AOM is insufficient to prevent methane from escaping to the bottom water, probably because of the abundant organic matter available for sulfate reducers to use instead of methane. In our sensitivity tests (using core 1CEL_Oct01 as an example), if AOM is allowed to be responsible for sulfate and methane consumption (no heterotrophic sulfate reduction and methanogenesis; $R_{\text{SD}}=R_{\text{AOM}})$ then methane concentrations would decrease to negative values (gray solid lines in Fig. 5b), which is inconsistent with observations. Although based on our data it is not possible to accurately quantify methane oxidation by calculating the relative proportion of sulfate loss due to heterotrophic sulfate reduction and/or AOM, our model results suggest AOM plays a minor role in this setting. It is also possible to have high rates of methane production and also AOM in the sediments but this is not captured as methane loss because there is more production than depletion. In addition the formation of gas bubbles and short residence time in the sediment due to shallow formation depth, low gas dissolution and fast release all contribute to lowering the relative impact of AOM. Other studies such as Lee et al. (2008) also detected the co-existence of porewater sulfate and methane in dwarf red mangrove habitats (Twin Cays, Belize) and and in that setting, like our study site, it was not possible to spatially separate the methanogenesis and AOM zones. Future investigations on the role of AOM in these dynamic mangrove-dominated tropical coastal lagoons are needed (e.g., Thalasso et al., 1997; Raghoebarsing et al., 2006; Lee et al., 2008; Kristensen et al., 2008; Beal et al., 2009; Silvan et al., 2011; Segarra et al., 2013).

Conclusions

The variable trends observed in sediment methane porewater chemistry from mangrove-dominated tropical coastal lagoons in Yucatán, Mexico, indicate a very dynamic system spatially and temporally throughout the year. This can be explained by multiple controlling parameters including physical processes such as mixing and dilution with seawater or groundwater, gas bubble rise and dissolution and microbial processes which operate at different rates during different times at all sites. Although our modeling suggests that organic carbon degradation rates are dominated by heterotrophic sulfate reduction in these cores, methanogenesis both in shallow and deeper sediments is also prevalent. The co-occurrence of methane and sulfate reduction (documented by sulfate depletion) in shallow sediments in this system is explained by high methane production rates supported by some combination of non-competitive substrates and ample dissolved and labile organic matter in the shallow sediments as well as the additions of methane from deeper sediment through gas rise and dissolution. Model results demonstrate that the largest sink for the methane in these sediments is efflux to the water column. Build-up of methane at shallow depths may reduce the fraction of methane that is oxidized prior to entering the water column, thereby increasing the flux at the sediment-water interface. This shallow methane pool may also encourage methane flux through bubble release, which can result in a larger fraction of the methane reaching the atmosphere without being lost to oxidation. Specifically, the ability of the microbial community in these sediments to use non-competitive substrates may allow for methane production in the upper sections of the sediment, potentially contributing to the higher than expected atmospheric methane flux measured from mangrove-dominated tropical coastal lagoons.

Modeling procedure used in the evaluation of porewater observations from sediments in mangrove-dominated tropical coastal lagoons, Yucatán, Mexico

Details of the modeling procedure and parameters used are described here. The following reactions are considered in the model:

Heterotrophic sulfate reduction ($R_{\text{SR}})$: $$\mathrm{2}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to \mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$$ Methanogenesis ($R_{\text{M}})$: $$\mathrm{2}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{CH}}_{\mathrm{4}}$$ Gas bubble dissolution ($R_{\text{MB}})$: $${\mathrm{CH}}_{\mathrm{4}}\left(\text{g}\right)\to {\mathrm{CH}}_{\mathrm{4}}\left(\mathrm{aq}\right)$$ The net reaction terms ($R_C$ in Eq. 2) are given in Table A1, boundary conditions are listed in Table A2, best-fit model parameters are given in Table A3 and model-derived concentration profiles are shown in Figs. 2 and A1.

 

[Figure omitted. See PDF]

 

[Figure omitted. See PDF]

Depth profiles of modeled (lines), measured (circles) and calculated (triangles) concentration of dissolved methane (dashed lines; open circles), sulfate (solid lines; solid circles) in the upper panel and sulfate depletion (solid lines; solid triangles), zero sulfate depletion (dashed lines) and chloride (gray circles) in the lower panel for each profile type (Groups 1–4, see text). One selected profile per group is shown in Fig. 2 for illustration and here the other profiles are shown (9 cores for Group-1, 6 cores for Group-2, 2 cores for Group-3 and 3 cores for Group-4). CEL and CH represent cores collected from Celestún Lagoon and Chelem Lagoon.

[Figure omitted. See PDF]

In Eq. (2), sediment porosity decreases with depth due to steady-state compaction: $$\mathrm{\Phi }={\mathrm{\Phi }}_f+\left({\mathrm{\Phi }}_{\mathrm{0}},-,{\mathrm{\Phi }}_f\right)\cdot e^{-px},$$ where ${\mathrm{\Phi }}_f$ is the porosity below the depth of compaction (0.78 for Celestún and 0.83 for Chelem), ${\mathrm{\Phi }}_{\mathrm{0}}$ is porosity at the sediment surface (0.90 for Celestún and 0.89 for Chelem) and $p$ ($\mathrm{1}/15$ cm${}^{-\mathrm{1}})$ is the depth attenuation coefficient. These parameters were determined from the measured porosity data at each site or at a nearby site (Eagle, 2002).

Under the assumption of steady-state compaction, the burial of porewater was calculated as in Berner (1980): $$v={\displaystyle \frac{{\mathrm{\Phi }}_f\cdot w_f}{\mathrm{\Phi }}},$$ where $w_f$ is the sedimentation rate of compacted sediments calculated from excess ${}^{210}$Pb data (0.25 cm yr${}^{-\mathrm{1}}$ for Celestún and 0.35 cm yr${}^{-\mathrm{1}}$ for Chelem; Gonneea et al., 2004). Sediment burial results in the downward movement of both sediment particles and porewater relative to the sediment water interface.

Rate expressions applied in the differential equations ($R_C$ in Eq. 2).

Variable	Rates	Applied cores
SO${}_{\mathrm{4}}^{\mathrm{2}-}$	$-R_{\text{SR}}$	Group-1
CH${}_{\mathrm{4}}$	$+R_{\text{M}}+R_{\text{MB}}$	Group-1
SO${}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}$	$+R_{\text{SD}}$	Group-1
SO${}_{\mathrm{4}}^{\mathrm{2}-}$	$+R_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$	Group-2, Group-3 and Group-4
CH${}_{\mathrm{4}}$	$+R_{{\text{CH}}_{\mathrm{4}}}$	Group-2, Group-3 and Group-4

Boundary conditions used in the model.

	SO${}_{\mathrm{4}}^{\mathrm{2}-}$	CH${}_{\mathrm{4}}$	SO${}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}$	SO${}_{\mathrm{4}}^{\mathrm{2}-}$	CH${}_{\mathrm{4}}$	SO${}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}$	Unit
	(top)	(top)	(top)	(bottom)	(bottom)	(bottom)
Group-1
1CEL_Apr00	5	0	4.8	8.5	0.5	5.534	mM
1CEL_Dec00	15	0.16	$-$2.2	5	0.56	2	mM
1CEL_Oct01	15	0	$-$2.3	7.5	0.295	4.6	mM
1CEL_Jul02	15	0.1	2.5	7.8	0.35	5.368	mM
2CEL_Jul02	18	0.02	10${}^{-\mathrm{9}}$	18.5	0.035	$-$2	mM
3CEL_Apr00	6.5	0.25	6.7	3.5	0.825	5.766	mM
3CEL_Jul02	13.8	0.31	2	6.5	1.3	3.5	mM
1_1CH_Oct01	15.1	0	12.4	13.2	0.0295	12.03	mM
1_2CH_Oct01	12	0.01	16	10	1	14.641	mM
2CEL_Dec00	21	0.01	$-$6.4451	7.6	0.25	4.6	mM
Group-2
1CH_Dec00	11.5	0.102		9.2	0.522		mM
1CH_Apr00	32	0.005		12.5	0.006		mM
2CH_Dec00	19.9	0.0015		7.96	0.0019		mM
5CH_Apr00	31.7	0.0031		29.1	0.0145		mM
2CEL_Oct01	5.0	0.511		7.88	0.734		mM
14CEL_Jul02	18.3	0.085		31.5	0.02		mM
16CEL_Dec00	8.8	0.038		8.81	0.025		mM
Group-3
5CEL_Apr00	17	0.047		11.6	0.0275		mM
14CEL_Dec00	20.5	2.1		34.9	0		mM
14CEL_Oct01	20	0.01		33	0.012		mM
Group-4
16CEL_Jul02	21	0		25.65	0.070		mM
16CEL_Oct01	23	0.00139		25.8	0.0015		mM
7CH_Oct01	20.5	0.00477		19.1	0.01		mM
8CH_Dec00	18.6	0		19	0.013		mM

Imposed and best-fit model parameters in each core.

	$T$	$S$	$P$	$D_{\mathrm{m}\left({{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}\right)}$	$D_{\mathrm{m}\left({\mathrm{CH}}_{\mathrm{4}}\right)}$	$D_{\mathrm{m}\left({{\mathrm{SO}}_{\mathrm{4}}}_{\mathrm{dep}}^{\mathrm{2}-}\right)}$	$L_{\text{MB}}$	$k_{\text{MB}}$	$k_{\text{SD}}$	$k_{{\text{CH}}_{\mathrm{4}}}$	$k_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$
	(${}^{\circ }$C)	(–)	(bar)	(cm${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}})$	(cm${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}})$	(cm${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}})$	(mM)	(yr${}^{-\mathrm{1}})$	(yr${}^{-\mathrm{1}})$	(yr${}^{-\mathrm{1}})$	(yr${}^{-\mathrm{1}})$
Group-1
1CEL_Apr00	27.3	17.6	1.06	354	598	354	1.2	1	500
1CEL_Dec00	22.2	16.4	1.06	367	523	367	1.3	0	400
1CEL_Oct01	31.2	13.9	1.1	382	659	382	1.4	0.6	500
1CEL_Jul02	30	21.1	1.01	374	640	374	1.1	0	500
2CEL_Dec00	22	17.7	1.06	315	520	315	1.3	1.6	500
2CEL_Jul02	28.7	20.8	1.01	364	619	364	1.1	0.1	500
3CEL_Apr00	28.6	20.2	1.07	363	618	363	1.2	0.9	400
3CEL_Jul02	30.4	18.2	1.01	377	646	377	1.1	50	500
1_1CH_Oct01	29.8	32.1	1.01	372	636	372	1.1	0	500
1_2CH_Oct01	29.8	32.1	1.01	372	636	372	1.1	0	500
Group-2
1CH_Dec00	25.2	24.8	1.05	318	556					1000	300
1CH_Apr00	26.3	39.4	1.09	347	583					1000	500
2CH_Dec00	23.9	27.5	1.08	329	547					4000	2000
5CH_Apr00	29.6	38	1.04	382	659					1000	1000
2CEL_Oct01	31.2	14.3	1.1	382	659					1000	500
14CEL_Jul02	31.5	27.4	1.01	385	663					1000	300
16CEL_Dec00	22.6	31.2	1.02	319	529					1000	300
Group-3
5CEL_Apr00	26.5	21.1	1.06	348	586					1100	3000
14CEL_Dec00	23.5	31.1	1.06	326	541					1000	300
14CEL_Oct01	31.1	13.9	1.07	382	657					1000	500
Group-4
16CEL_Jul02	30.3	30.5	1.01	376	644					600	300
16CEL_Oct01	29.7	28.2	1.06	319	529					1000	300
7CH_Oct01	29.6	31.3	1.01	371	633					500	300
8CH_Dec00	24.4	31.3	1.05	333	555					500	300

The sediment diffusion coefficient of each solute ($D_{\text{s}}$) was calculated according to Archie's law considering the effect of tortuosity on diffusion (Boudreau, 1997): $$D_{\text{s}}={\mathrm{\Phi }}^{\mathrm{2}}\cdot D_{\text{M}},$$ where $D_{\text{M}}$ is the molecular diffusion coefficient at the in situ temperature, salinity and pressure (Table A1) calculated according to Boudreau (1997). We used the same tortuosity coefficient (${\mathrm{\Phi }}^{\mathrm{2}}$ corresponding to $m=$ 3 in Archie's law) as reported by Wallmann et al. (2006) for fine-grained sediments.

Since net sulfate consumption is observed in Group-1 profiles (Figs. 2, A1), we used the following calculations to obtain net sulfate depletion rates ($R_{\text{SD}}$; mmol SO${}_{\mathrm{4}}^{\mathrm{2}-}$ cm${}^{-\mathrm{3}}$ yr${}^{-\mathrm{1}})$. $R_{\text{SD}}$ is proportional to the difference between modeled ($C\left({\text{SO}}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}\right))$ and measured concentrations ($C\left({\text{SO}}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}{)}_{\text{OBS}}\right)$: $$R_{\text{SD}}=k_{\text{SD}}\cdot \left(C,{\left({\mathrm{SO}}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}\right)}_{\text{OBS}},-,C,\left({\mathrm{SO}}_{{\mathrm{4}}_{\text{dep}}}^{\mathrm{2}-}\right)\right)$$ The corresponding kinetic constant is set to be high ($k_{\text{SD}}\ge$ 100 yr${}^{-\mathrm{1}})$ to ensure that simulated concentrations are very close to measured values. $R_{\text{SD}}$ implicitly includes $R_{\text{SR}}$ as well as anaerobic oxidation of methane ($R_{\text{AOM}})$: $${\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to {\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{HS}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$$ The numerical modeling procedure outlined in Wallmann et al. (2006) is used as a basis to simulate the rate of sedimentary organic carbon degradation ($R_{\text{POC}})$ by sulfate reduction and methanogenesis. Since the measured organic matter content in both lagoons showed evidence for a change in depositional pattern over time (Gonneea et al., 2004; Eagle, 2002), these measurements cannot be used for reliable organic matter degradation calculations. Hence, $R_{\text{SR}}$ (Eq. A5 below) was first calculated and then used to estimate $R_{\text{POC}}$ (Eq. A6) and subsequently to derive $R_{\text{M}}$ (Eq. A7). Here, we assume the three Reactions (R1), (R)2 and (R4) co-occur in the sulfate reduction zone such that the net reaction for methanogenesis and AOM (Reactions R2 and R4) is equal to carbon respiration by heterotrophic sulfate reduction (Reaction R1). In other words, $R_{\text{SD}}=$ 0.5$R_{\text{POC}}$.

To approximate the fraction of $R_{\text{POC}}$ due to $R_{\text{M}}$ and $R_{\text{SR}}$, a Michaelis–Menten kinetic limitation term is applied to Eqs. (A5)–(A7) (Wallmann et al., 2006): $$\begin{array}{}& & R_{\text{SR}}=R_{\text{SD}}=0.5\cdot R_{\text{POC}}\cdot f_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}& & R_{\text{POC}}={\displaystyle \frac{R_{\text{SD}}}{0.5\cdot f_{\mathrm{C}}\cdot f_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}}}& & R_{\text{M}}=0.5\cdot R_{\text{POC}}\cdot \left(\mathrm{1},-,f_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}\right)\end{array}$$ where $f_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}=\frac{C_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}}{C_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}+K_{\text{SR}}}$ is the Michaelis–Menten rate-limiting term for sulfate reduction.

At sites where methanogenesis was insufficient to simulate the measured methane data, methane was added as an external source by dissolution of gas bubbles (Chuang et al., 2013). Gas bubbles were observed in the field. The rate of dissolution of the gas bubbles (Reaction R3) rising through the sediment (CH${}_{\mathrm{4}\left(g\right)}\to$ CH${}_{\mathrm{4}})$ was also considered as (Haeckel et al., 2004) $$R_{\text{MB}}=k_{\text{MB}}\cdot \left(L_{\text{MB}},-,C_{{\text{CH}}_{\mathrm{4}}}\right)\text{if}{\mathrm{CH}}_{\mathrm{4}}\le L_{\text{MB}},$$ where $L_{\text{MB}}$ is the in situ methane gas solubility concentration calculated using the algorithm of Duan et al. (1992a, b) and the site-specific salinity, temperature and pressure. $R_{\text{MB}}$ depends on the first-order rate constant $k_{\text{MB}}$, which is a fitting parameter that lumps together gas dissolution in addition to diffusion of dissolved gas in the bubble tubes and walls.

Since sulfate depletion profile trends in Group-2, Group-3 and Group-4 show evidence of groundwater or seawater input with no positive depth integrated net sulfate depletion rates, the second approach for determining net methane and sulfate reaction rates for porewater data in these three groups is summarized as $$\begin{array}{}& & R_{{\text{CH}}_{\mathrm{4}}}=k_{{\text{CH}}_{\mathrm{4}}}\cdot \left(C_{{{\text{CH}}_{\mathrm{4}}}_{\text{OBS}}},-,C_{{\text{CH}}_{\mathrm{4}}}\right),& & R_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}=k_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}\cdot \left(C_{{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}_{\text{OBS}}},-,C_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}\right)\end{array}$$ Net methane and sulfate reaction rates are set to be proportional to the difference between modeled ($C_{{\text{CH}}_{\mathrm{4}}}$ and $C_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}})$ and measured concentrations ($C_{{{\text{CH}}_{\mathrm{4}}}_{\text{OBS}}}$ and $C_{{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}_{\text{OBS}}})$. The corresponding kinetic constants $k_{{\text{CH}}_{\mathrm{4}}}$ and $k_{{\text{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$ are listed in Table A3.

Methane fluxes at the boundaries were calculated using the model as follows: $$F_{{\text{CH}}_{\mathrm{4}}}\left(x\right)=\mathrm{\Phi }\left(x\right)\cdot \left(v,\left(,x,\right),\cdot ,C_{{\text{CH}}_{\mathrm{4}}},-,D_{\text{s}},\cdot ,{\displaystyle \frac{\text{d}{C}_{{\text{CH}}_{\mathrm{4}}}\left(x\right)}{\text{d}x}}\right),$$ where $x$ $=$ 20 cm is the bottom of the simulated core and $x$ $=$ 0 cm is the sediment–water interface.

Fixed concentrations were imposed for all solutes at the upper and lower boundaries to values measured at or near the sediment–water interface and at 20 cm. The method-of-lines was used to transfer the set of finite difference equations of the spatial derivatives of the coupled partial differential equations to the ordinary differential equation solver (NDSolve) in MATHEMATICA v. 7.0, using a grid spacing which increased from ca. 0.015 cm at the sediment surface to 0.38 cm at depth. Since most of the porewater profiles were fitted directly, only a few years of simulation time (5 years) was needed to achieve steady state. Mass balance was typically better than 99.9 %.

Acknowledgements

We thank the staff of the DUMAC Celestún station and the students of CINVESTAV for assistance with laboratory space, lodging, and field work. We thank Tom Lorenson and Ron Oremland of the Menlo Park, CA USGS for facility use and analyses for slurry incubations. This work was funded by Consejo Nacional de Ciencia y Tecnología Ref: 4147-P T9608, 32356T, and CONABIO Ref: B019 to Jorge A. Herrera-Silveira, NSF INT 009214214 to Adina Paytan, a Stanford Graduate Fellowship and Lieberman Fellowship to Megan B. Young and a fellowship of the Postdoctoral Research Abroad Program, sponsored by the National Science Council, Taiwan to Pei-Chuan Chuang (now: MOST; the Ministry of Science and Technology). We thank John Pohlman and an anonymous reviewer for their thoughtful comments and the associate editor (Helge Niemann) for handling the manuscript. Edited by: H. Niemann