File last modified 13 November 2000
18.4 An Example of CH4 at the FOAM Site
In this final section we will examine an example of 1-D modeling in anoxic sediments found in the Long Island Sound. This example is taken from a paper by Martens and Berner (1977, Limnol. Oceanogr., 22(1): 10-25). Of interest to these investigators, among many things, was the presence or lack of presence of SO4-2 where CH4 could be found.
The FOAM (Friends Of Anoxic Mud) site is located near-shore of Connecticut where the interstitial waters had incomplete to complete sulfate reduction in the upper ten or so centimeters. The overlying water depths were less than ten meters. The sediment cores were collected by gravity cores and the pore fluids were squeezed out by an inert gas-filtration system. The pore waters were analyzed for gasses and sulfate.
Of particular interest to these investigators was the small overlap of CH4 and SO4-2 in the muds they collected. At that time evidence had accumulated suggesting that methane did not form in sediments until all sulfate had been removed. In their sediment cores they found a small region of overlap between these two constituents (see Fig. 18.4.1) and wished to have one of four possible explanations tested. In the absence of evidence for actual bubble migration of methane into the overlying sediments Martens and Berner formulated the following hypothesis:
Martens and Berner noted that CH4 and SO4-2 did not seem to coexist together in these anoxic sediments, but there was some overlap implying that diffusion was responsible (CH4 from below and SO4-2 from above). The following reaction was suggested to account for the apparent mutual consumption:
which is thermodynamically favorable under the conditions at this location.
Given the lack of bioturbation in the sediment core and the a priori knowledge of the sedimentation rate (0.3 cm/yr) and methane diffusion coefficient (2 x 10-6 cm2/s) in this area from other work, they constructed a diagenetic equation to describe methane in the sediment column.
where k1 is the first order consumption rate (0.24 yr-1) methane estimated from "jar" experiments performed by other investigators and consistent with their own incubations.
While this equation has analytical solution, it is also an excellent place to introduce a numerical technique known as the tridiagonal algorithm.
The methane profile (above) is assumed to be in steady state and the concentration of methane is
"known" at the top and bottom of the profile. The choices for D, 
and k1 will determine the
shape of the profile between the end points (boundary conditions). This is a good example of the
two boundary point problem, but first the tridiagonal algorithm.
All differential equation's spatial differences can be expressed in interior-point form:
Let's take as an example the following generic equation:
which we can rewrite in FTCS form (without the advection divergence term):
Rearranging and multiplying through yields:
from which we can extract the following interior point coefficents:
As a general rule, round-off error is kept under control if Am>0, Bm>0, Cm>0 and Bm>Am+Cm. Now we can postulate the existence of two more vectors E and F such that:
and
Substituting 18.4.9 into 18.4.3 gives you:
Implying:
As you can see from these last two equations, all of Em and Fm can be expressed as a function of either the interior point coefficients (which we already know) and their own preceding values (that is to say, their m-1 points). From the lefthand side of 18.4.3 we can get the boundary conditions for E1 and F1 evaluate E and F recursively up to m=M-1. From here we can calculate Wm using 18.4.10 because at step m-1 Wm is known from the right hand side of 18.4.3's boundary condition and A, B, C, and D as well as E and F are already known or calculated.
An example should help clarify what is intended by this method.
For the methane profile shown in Fig. 18.4.1 we have the following equation (18.4.2 shown again):
We will apply the tridiagonal algorithm to this equation using the FUDM finite difference rendition (keep in mind that the tridiagonal algorithm can be applied to any finite difference form, FUDM is chosen because it provides a clear exposition).
Which we rewrite in another form as:
and since this problem is assumed to be at steady state, the lefthand side drops to zero. Now by
combining and rearranging into the interior point form we arrive at:
multiplying through by a negative one gives us our familiar interior point form and allows us to extract the following coefficients:
all of which we can evaluate. The trick now is to evaluate the boundary conditions and evaluate E1 , F1 and Wm and then put the machinery in motion. At the top of the sediment column (the seawater-sediment interface) the methane concentration is zero. Looking at 18.4.8 and 18.4.9 we can see that implies that E1=F1=0. At the other end of the profile we postulate that the methane concentration is at saturation (approximately 2600 micromolar), giving us Wm=2600.
This problem has been coded up in MatLab and is available in ch4_prof.m. You will need two other files, ch4.dat and tridiag.m, to run this program and see how it works. We encourage you to download it now and try it out.
GoTo Next Lecture