Biogeochemical consequences of macrofauna burrow ventilation†

The burrow walls created by macrofauna in aquatic sediments are sites of intense chemical mass transfer. Quantitative measurement of their significance is, however, difficult because chemistry in the immediate vicinity of burrow walls is temporally dynamic due to periodic ventilation of burrows by macrofauna. A temporally dynamic, 2D multicomponent diffusion-reaction model was utilized to depict the magnitude and time dependency of chemical mass transfer in the immediate vicinity of burrow walls as well as at the water/sediment interface. The simulation results illustrate that sediment particles, pore water, and microorganisms within a few millimeters of burrow walls experience significant oscillation in pH (as much as two pH units) and dissolved oxygen concentration (between saturation and near anoxia) whereas such oscillation is absent at the water/ sediment interface. The geochemical oscillation is expected to affect the net stability of mineral phases, activities and community structures of microorganisms, and rates and magnitudes of microbial diagenetic reactions.


Introduction
Burrowing infauna in aquatic sediments induce temporal fluctuation in concentrations of dissolved species (i.e., geochemical oscillation) through their metabolism and burrow ventilation activities. They periodically irrigate their burrows to replace metabolite-rich burrow water with fresh overlying water. The immediate vicinity of burrow walls is subject to periodic changes in the concentrations of oxygen, nutrients, and other pore water species. 1,2 Geochemical oscillation affects the courses of organic carbon (OC) diagenesis. 3 Field and laboratory studies suggest that OC remineralization rates are enhanced by the redox oscillation, even though the anoxic period in each oscillation cycle is typically much longer (y10-100X) than the oxic period for a given sediment microenvironment. 4,5 Temporally averaged redox conditions do not directly correlate to the rates and magnitudes of OC reactions. 3,6 A comprehensive understanding of OC diagenesis in aquatic sediments thus requires an adequate characterization of the temporal dynamics of redox and other geochemical parameters.
Directly observed data for geochemical oscillation induced by infauna are not abundant. A few studies have measured temporal dynamics of geochemical variables in biologically reworked sediments using redox potential, Eh, and oxygen microprobes. 1,7,8 The number of such studies is small despite the importance of geochemical oscillation: this may be largely due to the technical difficulties in measuring the geochemical phenomena that are not only temporally variable but also are confined to small regions. For example, in shallow marine and estuarine sediments, O 2 penetration in the vicinity of burrow walls is limited to a few millimeters, 9 and thus oscillation of redox and other related parameters may be restricted to the regions within a few millimeters of burrow walls. The temporal scale of oscillation at burrow walls depends on the burrow ventilation habits of infauna, and each oscillation cycle may be as short as several minutes. 10 Sampling devices typically used by geochemists (e.g., cores and pore water peepers) are designed to collect data that are averaged over the spatial and temporal extents too large to capture such variability. Microelectrode arrays and gel probes have elucidated sub-millimeter spatial variability of geochemical parameters, 11-13 but they have not been extensively used to capture the temporal variability.
Although the direct observations of geochemical oscillation in burrowed sediments are still few, such data can be put into the context of overall diagenesis when assisted with modeling. Recent studies have shown that the computational simulation of reaction couplings and kinetics associated with OC remineralization is a viable technique for the study of OC diagenesis. [14][15][16] Marinelli and Boudreau 11 applied this type of model to calculate the temporal fluctuation of O 2 and pH in the vicinity of model burrows artificially flushed with overlying water, and to estimate the effect of such oscillation on net chemical mass transfer.
In the present study, a numerical model of OC diagenesis and solute transport is used to investigate the magnitude of geochemical oscillation in the vicinity of burrows resulting from macrofauna ventilation and metabolite excretion. The model provides a tool for quantifying the variability of geochemical parameters that change rapidly within a small spatial extent.

Model geometry
The 2D cylinder geometry used in this study is taken from the single-component reaction, diffusion and burrow irrigation model first introduced by Aller. 17 The original model has since been expanded to a multi-solute numerical model for reaction, diffusion and discontinuous burrow irrigation. 11,18 The model geometry has also been adapted to allow burrows with variable depths. 19 The 2D cylinder geometry is advantageous to 1D vertical expression in describing temporal and spatial variability of geochemical parameters because it explicitly considers the radial chemical mass transfer in the vicinity of burrow walls. Whereas well-constrained 1D steady-state models are able to accurately hindcast net chemical mass transfer with the use of adjustable parameters, 20 they cannot represent radial variability such as the spatial and temporal variability seen in the immediate vicinity of burrow walls.
The model represents burrowed sediment with a series of equal cylindrical microenvironments in a closest packed array (Fig. 1). The outer radius of each cylinder (r 2 ) is the half distance between two neighboring burrows and is determined by the number of burrows per unit area of seabed. The model cylinder height is L, and the depth extent of all burrows are also assumed to be L. The concentration of a given solute species at a given spatial position (x, r) at time t, C x,r,t , can be determined by the following conservation equation that relates the time dependency of C x,r,t to vertical diffusion perpendicular to the water/sediment interface (WSI), radial diffusion perpendicular to burrow wall, and production or consumption of the solute due to biogeochemical reactions: 17,21 where Qwporosity, Dwdiffusion coefficient of the solute, hwdiffusive tortuosity, and Rwnet reaction rate. The diffusive tortuosity term in the equation can be replaced by a porosity expression using the following empirical correlation. 21 Diffusive tortuosity within the burrow cavity may be set to one, as the porosity therein is unity. The rate of OC remineralization (i.e., primary redox reactions, see Table 1) within burrow water may be set to zero since the amount of bacteria within the water is negligible compared to the amount within the sediments. However, for some solute species, metabolite excretion by the burrow inhabitants may be a significant contribution to the net chemical mass transfer. 22 Secondary redox reactions and acidbase reactions occur in both burrow water and sediment pore water.

Boundary conditions
Boundary conditions needed for the numerical solution of eqn. (1) include the following: 11 Solute concentrations at WSI (i.e., x~0) are assumed to be constant (C 0 ) (eqn. (3)). Solute concentrations of the burrow water during the ventilation period are also assumed equal to the overlying water (eqn. (4)), whereas, during the rest period, the burrow water composition is determined according to eqn. (1) by diffusive transport across the burrow wall, macrofaunal metabolite production, and aerobic reoxidation of reduced species. In reality, burrow water composition does not remain identical to that of the overlying water during ventilation because the volume of water replaced by each pumping of macrofauna is smaller than the volume of burrow cavity and the velocity of water due to macrofauna flushing has a finite value. 22,23 Radial solute flux at the outer surface of each cylinder (i.e., r~r 2 ) and at burrow axis (i.e., r~0) are set to zero due to symmetry (eqn. (5) and (6)). Vertical solute flux at the bottom boundary (at x~L) is set to zero (eqn. (7)). The model geometry is based on the assumption that all burrows have the depth extent of L. In actual simulations, L has to be set greater than the depth extent of one's interest in order to avoid numerical artifacts. This is obviously unrealistic because actual burrows have variable vertical extents and some of them may be much shorter than L. Consequently, the model is expected to be the most accurate near WSI: its accuracy decreases with depth. At the burrow wall (i.e., r~r 1 ), solute concentrations and radial flux are both continuous (eqn. (8) and (9)).

Solution scheme
A FORTRAN code is written to numerically solve the conservation equation (1) with both left and right hand sides of the equation being fully discretized. A 2D uneven grid 24 is utilized in order to have small grid spacings near WSI and burrow wall (Dx~Dr~0.25610 23 m) where rapid aerobic reactions are expected to yield steep spatial gradients in concentrations. The C x,r for each grid point is solved as a time-evolution problem until the evolution of C x,r distribution during a given ventilation-rest cycle becomes identical to the evolution during the previous ventilation-rest cycle. The conservation equation is written for each of the following species: O 2 , NO 3 2 , SO 4 22 , NH 4 z , SS (wH 2 SzHS 2 ), SCO 2 (wCO 2 zH 2 CO 3 zHCO 3 2 zCO 3 22 ), and titration alkalinity (Alk t~H CO 3 2 z2CO 3 22 zHS 2 ). All conservation equations are coupled and solved simultaneously at each time step through the reaction terms, R, as shown in Table 1. For example, the reaction term for the conservation equation of O 2 at time step T is determined by the concentration values of O 2 , NH 4 z , and SS at time step T21 (see eqn. (I-1)). The diffusion coefficients of solute species, Monod kinetics constants, and thermodynamic constants for sulfide and carbonate systems can be taken from previously established equations found in literature that express these parameters as functions of temperature and salinity. 21 The other siteand species-specific parameters, including the parameters for reaction kinetics and burrow geometry, need to be taken from each environment being studied. The durations of ventilation and rest periods are also required to implement the periodic, discontinuous irrigation scenario.

Results
Directly determined data on macrofauna burrow geometry, ventilation habits, OC degradation rates, porosity, and pore water chemistry are necessary in order to properly constrain the model and evaluate the simulation results. Collective studies of sediments populated with Nereis diversicolor 25-28 provide many of the necessary data.

Model system description
Model simulations are carried out for the mesocosm systems of Kristensen and Hansen. 25 All necessary parameters were given in Kristensen and Hansen 25 or the references therein. [26][27][28] The mesocosms were loaded with homogenized fine-grained mud with 75% porosity to the depth of 15 cm, and populated with 1200 m 22 Nereis diversicolor of 200-400 mg wet weight. Each N. diversicolor constructed a U-shape burrow with two openings at WSI. The upper 5-6 cm of the sediment columns were burrowed. The ratio for the production of SCO 2 and NH 4 z due to microbial remineralization was given by Kristensen and Hansen 25 to be 4.8, which they determined by correlating the pore water SCO 2 and NH 4 z profiles using the method previously described. 28 These values are much lower than the C : N ratio of bulk organic matter in the source mud from Kertinge Nor, Denmark (#12) given by Hansen and Kristensen 28 because organic matter that is readily decomposed by the microbial remineralization has a lower C : N ratio than the bulk organic matter pool. 28 The bottom water SCO 2 and NH 4 z concentrations were measured to be 1.82610 23 M and 0.014610 23 M, respectively. Kristensen and Hansen 25 also described control mesocosms that were loaded with the same mud but with no macrofauna.
The anaerobic OC degradation rate, R OC an , was determined by fitting the 1D steady-state model to measured SCO 2 and NH 4 z depth profiles of control mesocosms. The control mesocosms were not burrowed, thus the model used for this simulation was 1D, rather than 2D, and assumed steady state: This 1D conservation equation was written for each of the pore water species shown in Table 1. The equations were then Table 1 Reactions and rate expressions considered in the model simulations. Rate symbols are defined in Table 2 Primary redox reactions- Secondary redox reactions- Reaction rates- coupled and solved numerically. The simulation determined the best-fit R OC an value to be 2.5610 29 M s 21 (Fig. 2). The R OC an value was assumed to be identical for both control and irrigated mesocosms because no data are available that would allow a reasonable estimate of rate enhancement due to irrigation. All other parameters were directly determined using actual data in the original study, or were taken from well-established literature data. The parameter values used are shown in Table 2. In these calculations, the surface adsorption of NH 4 z was not considered: this is because the C : N ratio of remineralization given by the original study 25 (4.8) was calculated without considering adsorption, thus including implicitly the effects of NH 4 z loss due to sorption. 28 R OC an was treated as depth independent because the experimental mud was homogenized prior to loading. 25 Using the parameter values shown in Table 2, the cylinder model was constrained. The additional parameters necessary to properly constrain the temporally dynamic 2D model were taken from the original studies of N. diversicolor. [25][26][27] The radius and volume of N. diversicolor burrows were represented by 1.0610 23 (m) and 0.554 (cm 3 ), respectively, by taking the median values of burrow diameters and lengths estimated for medium-sized N. diversicolor. 31 The durations of ventilation and rest were taken from the observations by Miron and Kristensen. 26 The NH 4 z excretion rate of N. diversicolor was taken from the study of Nithart et al. 27 in which the rate was measured in an estuarine environment at 15 uC to be 0.0111 (mM h 21 mg 21 dry wt.). Using the dry/wet wt. ratio for average polychaetes 32 (~0.13), and the median wet wt. of N. diversicolor in the original study of Kristensen and Hansen 25 (300 mg) and median burrow volume (0.554 cm 3 ), the NH 4 z excretion rate of N. diversicolor into the burrow water of experimental mesocosms was estimated to be 2.18610 27 M s 21 . Using the same C : N ratio as the microbial remineralization of the mesocosm mud (4.8), the SCO 2 excretion to burrow water by N. diversicolor was consequently determined to be 1.05610 26 M s 21 . In reality, the C : N ratio of organic matter metabolized by macrofauna may be different from the C : N ratio of microbial remineralization for a given mud. The parameter values are summarized in Table 3.

2D simulation results
Two-dimensional model simulations depict the oscillation of O 2 and pH in the immediate vicinity of WSI and burrow walls during each N. diversicolor rest-ventilation cycle ( Fig. 3 and 4). The simulation results indicate that the sediment particles and microorganisms that are located at burrow walls experience oscillations in O 2 concentration (2.3610 24 to 1.1610 24 M) and pH (8.0 to 6.9) during each 556-s period. The burrow walls are marked by the steep gradients in geochemical parameters that include pH and O 2 . Whereas WSI also is a site of steep geochemical gradients, it is temporally stable because it interfaces the overlying water whose composition remains constant. The magnitude of pH decrease during the rest period is dictated by the SCO 2 increase in burrow water due to N. diversicolor metabolism (Fig. 5).
The 2D simulation results were also used to calculate radially averaged depth profiles of NH 4 z and SCO 2 . The outcome is compared to the profiles measured in the mesocosms studied by Kristensen and Hansen 25 (Fig. 6). Temporal variability in the profiles during each rest-ventilation cycle is insignificant when the concentration values are radially averaged, and thus only   The relief in these mesh diagrams indicates the concentrations and gradients. At the beginning of each rest period (T1~0 min), the burrow water is fully oxygenated. It becomes less oxygenated during the rest period, as seen at halfway into the rest period (T2~2.3 min) and at end of the rest period (T3~4.6 min). Immediately following, the macrofauna begins ventilation (T4~4.6zD min; D is equal to one time step in the model run, which is 8 in these simulations) resulting in the maintenance of the burrow water O 2 level, as seen at halfway into the ventilation period (T5~6.9 min) and at the end of each rest-ventilation cycle (T6~9.3 min). Fig. 4 Model-calculated evolution of pH near WSI and burrow wall during each rest-ventilation cycle of N. diversicolor. At the beginning of each rest period (T1~0 min), the pH along burrow wall is 8.0. It decreases during the rest period, as seen at halfway into the rest period (T2~2.3 min), and at the end of the rest period (T3~4.6 min) when the pH at burrow wall is as low as 6.9. Immediately following the rest period, the macrofauna begins ventilation (T4~4.6zD min; D is equal to one time step in the model run, which is 8 in these simulations) resulting in the maintenance of the burrow water pH level, as seen at halfway into the ventilation period (T5~6.9 min) and at the end of each rest-ventilation cycle (T6~9.3 min). one profile is shown for each species. The simulated depth profiles agree with the measured profiles for the upper 1-2 cm of the mesocosm sediments where the actual number of burrows is well represented by the model geometry. The agreements become much less satisfactory at depths where there are fewer burrows than represented by the model geometry, which assumed all N. diversicolor (1200 m 22 ) to construct vertical burrows of the uniform radius with two openings as deep as the height of the mesocosms. The disagreement may also be due to the assumption of 100% irrigation during the ventilation period. In reality, the composition of burrow water is not identical to that of overlying water during the ventilation period, and the difference between the burrow water and overlying water may be greater in the deeper part of the burrow than in the near-surface part of the burrow.

Effect of ventilation patterns
During the rest period, the burrow water becomes less oxygenated, and more enriched in metabolites such as NH 4 z and SCO 2 . The build up of SCO 2 also induces decrease in pH values. The above calculations imply that the duration of the rest period is an important variable in determining the magnitudes of geochemical oscillations. The longer the rest period, the lower the pH values and O 2 concentrations would become, and thus the ranges of oscillation during each restventilation cycle would become greater. In order to illustrate the significance of these variables, the 2D model was applied to a model system that consists of the same mud as the original mesocosms of Kristensen and Hansen 25 but is populated with 1200 m 22 N. virens instead of N. diversicolor which has the rest period of t rest~2 391 s and ventilation period of t vent~8 40 s, according to a mesocosm study by Kristensen and others. 22 This rest period is nine times greater than the rest period of N. diversicolor.
N. virens construct U-shaped burrows. The model burrow diameter (r 1 ) was determined from the mesocosm study to be r 1~3 .0610 23 m. 22 The same study 22  The results of N. virens simulation ( Fig. 7 and 8) indicate that the sediment particles and microorganisms that are located at burrow walls would experience significant oscillations in O 2 concentration (2.3610 24 to 0.1610 24 M) and pH (8.0 to 5.8) during each 3231-s rest-ventilation cycle. The magnitude of oscillation is greater in this N. virens scenario than for N. diversicolor because the longer rest period promotes build up of metabolites such as NH 4 z and SCO 2 due to macrofaunal metabolism (Fig. 9) and loss of O 2 due to diffusive loss and aerobic reoxidation of reduced species. Kristensen 2 observed Fig. 5 Model-calculated evolution of SCO 2 near WSI and burrow wall during each rest-ventilation cycle. During the rest cycle (T1-T3), SCO 2 within burrow cavity evolves to exceed that of adjacent sediment pore water due to metabolic SCO 2 excretion by N. diversicolor. complete anoxia within N. virens burrows within 5-7 min of the onset of rest period. On the other hand, the model burrow of this study does not reach anoxia during the rest period because of the assumption that the burrow water is fully oxygenated during the ventilation period. In reality, O 2 concentration of burrow water remains lower than that of overlying water even during the ventilation period because each pumping by macrofauna is likely to replace only part of the water in the Fig. 7 Model-calculated evolution of O 2 concentration near WSI and burrow wall during each rest-ventilation cycle of N. virens (54 min). At the beginning of the rest cycle (T1~0 min), the burrow water is fully oxygenated whereas it becomes depleted with O 2 during the rest period, as seen at halfway into the rest period (T2~20 min) and at the end of the rest period (T3~40 min). Immediately following the rest period, the macrofauna begins ventilation (T4~40zD min; D is equal to one time step in the model run, which is 8 in these simulations) resulting in the maintenance of the burrow water O 2 level, as seen at halfway into the ventilation period (T5~47 min) and at the end of each rest-ventilation cycle (T6~54 min). burrow cavity. Macrofauna itself may also consume O 2 still remaining in its burrow during the rest period.

Difference between two types of interface
The simulations indicate that both WSI and burrow walls are the sites of intense chemical mass transfer as evidenced by the steep gradients in geochemical parameters such as pH and O 2 concentrations. Burrowed sediments accommodate greater extent of aerobic OC remineralization than non-bioturbated sediments due to the increased interfacial area between anoxic sediments and oxygenated water. There is, however, one important difference between the redox interface at WSI and the redox interface at burrow walls. The burrow walls are temporally dynamic because they border burrow water whose geochemical properties oscillate according to the metabolism and ventilation of burrowing macrofauna. On the other hand, WSI is adjacent to the overlying water whose geochemical properties are temporally more stable. Consequently, sediment particles and microorganisms in the close vicinity of burrow walls experience temporal oscillation in redox chemistry and pH whereas sediment particles and microorganisms in the immediate vicinity of WSI do not undergo such oscillation. For example, microorganisms in the close vicinity of N. virens burrows have to adapt to the oscillating redox environment in which they may need to switch between different terminal electron acceptors for respiration during a short period of time. In addition, the same microorganisms are required to tolerate the oscillation in pH, which has been shown to influence the growth rates and activities of bacteria. 33 The mineral particles near the burrow walls also experience the same pH oscillation. If mineral particles such as calcareous tests of microorganisms or iron oxyhydroxides are located in the immediate vicinity of burrow walls, they may experience oscillation between supersaturation and undersaturation due to the fluctuation in pH.
These two different types of interfaces (WSI and burrow walls) are expected to accommodate different microbial ecology and distribution of functional groups because their redox, nutrients, and pH environments are distinctively different (i.e., steady state vs. periodic oscillation). While the enhancement of net microbial activities due to the presence of burrows is well documented, 34,35 the mechanisms for enhancement have not been fully evaluated. Potential enhancement mechanisms include the faster removal of toxic metabolites, increased diffusive transport, increased introduction of terminal electron acceptors, and resulting stimulation of production. 35 Further investigations of the spatial and temporal variability of microbial activities, community structures, and geochemical environments in the immediate vicinity of burrow walls are warranted.

Conclusions
Numerical simulations using a 2D diffusion-reaction model were conducted to examine the geochemical effect of macrofaunal burrow ventilation in the immediate vicinity of burrow walls and at the water/sediment interface. During the rest period, macrofaunal metabolites such as NH 4 z and SCO 2 build up while loss of O 2 occurs due to molecular diffusion and aerobic reoxidation of reduced species within the burrow cavity water. On the other hand, the geochemical composition of burrow cavity water becomes rich in O 2 and depleted in metabolites during the ventilation periods. Consequently, the redox and other geochemical parameters, including pH, oscillate in the immediate vicinity of burrow walls during each rest-ventilation cycle. The range of oscillation is greater when the burrow occupant maintains longer rest periods. Such oscillation is absent in the vicinity of the water/sediment interface, which borders the overlying water whose composition remains constant.
The redox and pH oscillations depicted through the model simulations are likely to influence the activity and community structure of sedimentary microorganisms and thermodynamic stability of mineral particles. Further studies are necessary in Fig. 9 Model-calculated evolution of SCO 2 concentration near WSI and burrow wall during each rest-ventilation cycle of N. virens. Significant increase in the metabolite concentrations occurs toward the end of rest cycle (T3~40 min). order to quantify the couplings between geochemical variables and microbial properties in temporally dynamic environments such as the vicinity of burrow walls.