Quantifying bioirrigation using ecological parameters: a stochastic approach†

Abstract

Irrigation by benthic macrofauna has a major influence on the biogeochemistry and microbial community structure of sediments. Existing quantitative models of bioirrigation rely primarily on chemical, rather than ecological, information and the depth-dependence of bioirrigation intensity is either imposed or constrained through a data fitting procedure. In this study, stochastic simulations of 3D burrow networks are used to calculate mean densities, volumes and wall surface areas of burrows, as well as their variabilities, as a function of sediment depth. Burrow networks of the following model organisms are considered: the polychaete worms Nereis diversicolor and Schizocardium sp., the shrimp Callianassa subterranea, the echiuran worm Maxmuelleria lankesteri, the fiddler crabs Uca minax, U. pugnax and U. pugilator, and the mud crabs Sesarma reticulatum and Eurytium limosum. Consortia of these model organisms are then used to predict burrow networks in a shallow water carbonate sediment at Dry Tortugas, FL, and in two intertidal saltmarsh sites at Sapelo Island, GA. Solute-specific nonlocal bioirrigation coefficients are calculated from the depth-dependent burrow surface areas and the radial diffusive length scale around the burrows. Bioirrigation coefficients for sulfate obtained from network simulations, with the diffusive length scales constrained by sulfate reduction rate profiles, agree with independent estimates of bioirrigation coefficients based on pore water chemistry. Bioirrigation coefficients for O₂ derived from the stochastic model, with the diffusion length scales constrained by O₂ microprofiles measured at the sediment/water interface, are larger than irrigation coefficients based on vertical pore water chemical profiles. This reflects, in part, the rapid attenuation with depth of the O₂ concentration within the burrows, which reduces the driving force for chemical transfer across the burrow walls. Correction for the depletion of O₂ in the burrows results in closer agreement between stochastically-derived and chemically-derived irrigation coefficient profiles.

Introduction

Nearshore sediments support large populations of burrowing macroinfauna, such as polychaete worms, shrimps and crabs. These organisms build extensive burrow networks within the upper 10 to 100 cm of the sediment column, which are either passively or actively flushed with water from the overlying sediment/water interface (SWI). Differences in solute concentrations between water flushed through burrows and pore water in the surrounding bulk sediment provide the driving force for enhanced solute mass transfer, called bioirrigation. [1–3] Bioirrigation not only enhances solute fluxes,[4, 5] but also influences the microbial community structure, [6–9] as well as sediment and pore water composition. [10–15]

Mathematical one-dimensional (1D) models of bioirrigation have been developed that treat bioirrigation variously as a diffusive, advective or nonlocal process. Diffusive models of bioirrigation use increased diffusion coefficients to account for enhanced solute transport due to macrofaunal activity within a defined mixing zone. [16–20] A diffusive representation implies that biological solute mixing is rapid and occurs on a spatial scale much smaller than the characteristic length scale of chemical pore water gradients.[21] Less common are advective models of bioirrigation, in which an enhanced advection velocity accounts for increased vertical pore water transport. [11, 22, 23] Boudreau[24] and Boudreau and Imboden[25] have argued that a more appropriate description of bioirrigation is based on a nonlocal transport description. The most common approach employs a mass transfer coefficient or bioirrigation coefficient, α, to describe the rate of exchange between overlying water and pore waters at depth. [12, 26, 27] None of the 1D bioirrigation models, however, include an explicit representation of the burrow networks constructed by the resident macrofaunal populations.

An important step in bridging the gap between chemical transport modeling of bioirrigation and benthic ecology was the development of the 3-dimensional (3D) radial diffusion model by Aller.[1] In this approach, burrows are represented as continuously flushed, evenly spaced cylindrical burrows of constant diameter and depth. Solutes are transported diffusively across the burrow walls in response to concentration gradients between bulk pore waters and waters flushed into burrows from the overlying SWI. Boudreau and Marinelli extended the 3D burrow model to describe discontinuous flushing of burrows, as well as variations in animal population (distance between burrows) and organism sizes (burrow radii).[5] Boudreau further showed that the 1D nonlocal description of irrigation is equivalent to the 3D model of continuously flushed vertical cylinders; the 1D nonlocal equation is simply the radially integrated form of the 3D model.[28] The equivalence between the 3D burrow model and the 1D nonlocal model implies that the magnitude of bioirrigation at a given depth depends critically on the burrow wall surface area at that depth.

In all of the above models (diffusive, advective, nonlocal and cylindrical burrow) assumptions are made regarding the depth dependence of the intensity of bioirrigation. Bioirrigation intensity is either assumed to be constant over a discrete interval, or a functional depth-dependence is assigned to enhanced diffusion coefficients or nonlocal exchange coefficients.[12, 13, 29, 30] Values of irrigation coefficients are generally derived by fitting chemical profiles, rather than by relying on ecological information. However, in principle, the latter could be used to determine the depth dependence of bioirrigation intensity directly. In addition, the models discussed so far are strictly deterministic and provide no information concerning the spatial or temporal variability of irrigation.

In this study, a stochastic model of bioirrigation based on benthic ecological data is presented. Information on burrow sizes and shapes for a variety of representative macrofaunal organisms is used to calculate stochastic realizations of 3D burrow networks. By generating a large number of burrow networks it then becomes possible to determine both average burrow surface areas and their variability as a function of depth. Benthic ecologists have characterized the structures of burrows using methods such as resin casting and X-radio-graphy; these detailed measurements can be used to calculate the depth-dependence of burrow surface areas without implementing a stochastic approach. However, the stochastic model presented here is advantageous because it can also be used to provide estimates of surface area variability. Surface areas from the stochastic model provide an independent method for constraining the depth dependence of irrigation coefficients used in 1D nonlocal models of bioirrigation. The stochastic model also provides a method to assess the uncertainty associated with irrigation coefficients in sediments.

Model

The stochastic burrow network

In the stochastic model, the sediment is represented as a 3D grid over which burrows are distributed. Burrows are modeled using ten endmember shapes commonly observed in macrofaunal burrow networks (Fig. 1). The following input parameters, with corresponding standard deviations and probability distributions where appropriate, are user-specified: (1) total burrow density (burrows per m²), (2) probability of occurrence for each of the 10 possible burrow shapes, and (3) radius and lengths of segments associated with each burrow shape. For example, U-shaped burrows require specifications for vertical and horizontal segment lengths; Y-shaped burrows require specifications for lengths of inclined upper branches and for the horizontal stem.

Figure 1

Endmember shapes used in stochastic burrow network simulation.

With the specified input parameters, burrows are distributed over the grid as follows. For each realization of the network a total burrow density is chosen randomly from the user-specified probability distribution of the total burrow density. The total burrow density is partitioned among the different burrow shapes based on the user-specified probabilities of occurrence of the ten possible burrow shapes. For each burrow, the lengths of the burrow segments and the radius of the burrow are selected similarly using the corresponding probability distributions. For nonvertical burrow orientations, the direction of the horizontal or inclined burrow section is selected randomly as north, south, east or west, with equal probability for each direction.

After all of the burrow parameters have been calculated, the burrow is placed into the 3D network of grid blocks. Each block either contains bulk sediment or is part of a burrow. The location of the burrow is determined by generating random coordinates for the burrow exit at the top surface of the grid (i.e., at the SWI). For each burrow type, the user specifies whether the burrow is allowed to intersect with other burrows. If the calculated exit coordinates cause the burrow to intersect with another burrow, and no intersection is allowed, the exit coordinates are recalculated until the burrow can be located in the grid, or until a large number of attempts (specified by the user) is exceeded, ending the simulation. To eliminate edge effects, burrows that exit the grid sides are continued at the opposite edge (i.e., a periodic boundary condition is imposed). For all stochastically generated parameters, either a Rayleigh or Gaussian probability distribution function may be selected, or the parameter may be specified deterministically for use in all simulations.

For each stochastic realization of the burrow network (Fig. 2), the burrow volume, density and burrow wall surface area are calculated within each depth interval of the 3D grid. To calculate surface areas, all burrows are treated as cylindrical tubes. The surface area (s_h) per unit sediment depth of vertical or horizontal burrows is therefore equal to

Figure 2

Sample 3D stochastic network simulation of burrows and corresponding burrow surface area per volume (red line) as a function of depth.

where r₁ is the burrow radius, β is the angle of the burrow with respect to the vertical, and the subscript h indicates that the surface area within the depth interval is normalized to the uniform grid spacing, h. The total burrow surface area per unit sediment depth, S_h, at any given depth is then obtained by:

S_h = ∑ s_h     (2)

where summation is carried out over all the burrows in the given depth interval. For simplicity, in the case of intersecting burrows, the surface area at the grid point of intersection is that calculated for the last burrow placed in the grid block. Because burrows are placed into the grid randomly and many networks are simulated, this results in an averaging of intersecting large and small surface area burrows. The mean values and standard deviations of the burrow density and burrow surface area are calculated as a function of depth by performing a large number of successive burrow network simulations (default 10 000 simulations).

Burrow networks and nonlocal exchange coefficients

Derivation of nonlocal exchange coefficients from burrow surface areas

Bioirrigation occurs because burrows provide conduits for solute mass transport that is rapid relative to diffusion through the bulk sediment. The burrow network therefore significantly influences biogeochemical cycling in the sediment. This is especially true in nearshore environments, where biologically mediated transport is often the dominant solute transport mechanism in surface sediments.[31] Studies that attempt to link observed pore water profiles with the underlying process dynamics thus require a mathematical description of solute exchange due to bioirrigation. One such model is the nonlocal exchange coefficient description of bioirrigation,[26] where the rate at which a solute is added or removed by bioirrigation (I_n1) at a given depth in the sediment is given by

I_nl(x) = α _i φ(C₀ - C_avg)| _x ,     (3)

where φ is the sediment porosity, and the bioirrigation coefficient, α _i , in units of inverse time, denotes the nonlocal exchange coefficient describing mass transfer of a given solute, i, between the overlying water and the bulk sediment; C₀ and C_avg are the solute concentrations in the overlying water and in the bulk sediment, respectively.

Boudreau has demonstrated that the nonlocal exchange coefficient description of bioirrigation is equivalent to the continuously flushed radial tube model of bioirrigation developed by Aller. [1, 28] He proposes the following equation:

where D _i is the molecular diffusion coefficient of the i-th species corrected for tortuosity, r₁ the radius of the vertical, cylindrical burrows, r₂ the half-distance between burrows, and (∂C/∂r) the concentration gradient of solute i extending from the burrow wall into the bulk sediment. In the idealized radial tube model, the term 2r₁/(r₂²-r₁²) corresponds to the burrow surface area per volume of sediment surrounding an individual burrow (s_v). The main difficulty in applying eqn. (4) resides in estimating (∂C/∂r). If we approximate the concentration gradient by a truncated Taylor series expansion, then

where

is the (time averaged) distance from the burrow wall at which the solute concentration reaches the radially averaged concentration, C_avg, for the bulk sediment at a given depth and C_b, is the (time averaged) concentration in the burrow.

If the total burrow volume is relatively small (i.e., r₁ ≪ r₂), the volume of bulk sediment surrounding an individual burrow at a given depth is approximately equal to the total sediment volume divided by the number of burrows at that depth. Therefore, s_v ≈ S_v, where S_v is the burrow surface area summed over all burrows divided by the total volume of sediment modeled within a given depth interval. At the sites assessed in this study (Dry Tortugas, FL and Sapelo Island, GA), the contribution of burrows to total sediment porosity is small (< 5%), justifying this simplification. If we assume that the solute concentrations within burrows are similar to the concentrations in the overlying water (C_b, ≈ C₀) then, using eqn. (3)-(5), it follows that

Eqn. (6) is derived for the idealized model of evenly spaced, identical vertical burrows. In this study, we assume that, within each horizontal slice of the 3D sediment grid, eqn. (6) provides a good approximation for the relationship between the bioirrigation intensity and the size and density of burrows. The stochastic burrow network can then be used to calculate mean values of burrow surface areas and radii, which can be implemented in eqn. (6). While the assumption underlying the equivalence between 1D nonlocal irrigation coefficients and idealized 3D cylindrical burrows is not strictly met for complex burrow geometries, the approach provides a means to estimate the depth-dependence of α, which is otherwise poorly constrained.

The equivalent burrow radius r₁(x) is calculated from the total burrow surface area per unit sediment depth according to

where n_tot(x) is the total number of burrows encountered in the depth interval over which S_h is determined. Both n_tot and S_h are calculated from the stochastic realizations of the 3D burrow network. Eqn. (7) accounts for the fact that inclination of the burrows with respect to the vertical increases the surface area available for exchange between burrows and bulk sediment.[31, 32] The average half-distance between burrows may be required for estimation of (see below). It is obtained by analogy with an ideal closest packing arrangement of vertical burrows, i.e. burrows (and the cylindrical packet of sediment surrounding each burrow) are assumed to be distributed evenly over the grid. The average half-distance between burrows at a given depth x, is then given by

where n_tot(x) is the mean total number of burrows encountered in the sediment slab, and A₀ is the horizontal surface area of the 3D grid.

Diffusion length scales

From eqn. (6) it is evident that the irrigation coefficient may depend on the chemical species, not only through the species-specific diffusion coefficient, but also through a characteristic radial diffusion length scale, L _i = - r₁. In theory, the latter can be obtained by direct measurements of chemical microprofiles across at the burrow/sediment interface (BSI). Such data, however, are rarely available. Alternatively, vertical microprofiles recorded at the sediment/water interface (SWI) may be used as analogs for the gradients at the BSI, or the radial diffusion length scales may be estimated from simplified reaction-transport models.

If vertical diffusion at the SWI is considered analogous to radial diffusion at the BSI, the measured vertical length scale of the concentration change at the SWI (

) must be corrected to convert from the planar geometry of the SWI to the radial geometry of the BSI. If we assume that the net rate of consumption (production) of the solute of interest is similar in the vicinity of the SWI and the BSI, then we can equate the planar and cylindrical volumes in which the solute is consumed (produced). Specifically, at the SWI, the ratio of a volume of a unit area of surface extending to depth , divided by the surface area of the sediment is . This must be equal to the ratio of sediment volume associated with an individual burrow to the burrow's surface area,

which yields

The radial diffusion length scale

predicted by eqn. (10) is strictly speaking only valid near the SWI, where reflects the local production or consumption of the solute. As the production (consumption) rate changes with depth, so will the value of .

When measurements of rates of consumption or production of the solute in the bulk sediment are known, it is possible to estimate values of L _i by assuming that in close proximity of burrows solute transport is dominated by irrigation (i.e., vertical diffusion can be ignored). From the vertical profiles of solute concentration and net rate of production/consumption α can be determined by balancing the flux across the burrow walls with the reaction taking place in the surrounding sediment,

Using the approximation of eqn. (5), and assuming C_b ≈ C₀, the corresponding length scale at the BSI for r₂ ≫ r₁ is then

where R _i is the net rate of production/consumption in the bulk sediment expressed per unit pore water volume.

Estimates of length scales and bioirrigation coefficients can also be obtained if the reaction kinetics are known and the governing mass balance equation is solved at steady state, with the average concentration defined as

Subsequently,

can be determined using the condition C() = C_avg. For example, for zeroth order kinetics,

With the boundary conditions C(r₁) = C₀ and (∂C/∂r) at r₂ = 0, the radial concentration profile is given by

Solving for

leads to

where f(x) fulfills f(x)·e^f(x) = x. For |x| < 0.3, f(x), the Lambert W function, can be expressed with a relative error <1% in a series expansion,

From eqn. (16) it is apparent that for zeroth order kinetics,

depends only on the burrow size and spacing, not on the reaction rate or diffusion coefficient. Thus, the expected species specificity of α only results from the species dependence of the effective diffusion coefficient. For other reaction kinetics, [1, 33] numerical methods are required to determine .

Burrow flushing efficiency

Nonlocal transport models of irrigation generally assume perfectly efficient flushing of burrows, that is C_b = C₀ (eqn. (3), (6). If perfect flushing does not occur, then depletion (or accumulation) of the solute within the burrow will cause the nonlocal irrigation coefficient, α_i,b(x) to deviate from the apparent irrigation coefficient, α_i,0(x) defined by eqn. (3). The relationship between the true irrigation coefficient (α_i,b) and the apparent coefficient is

For highly reactive species, or small burrows, this effect may significantly affect bioirrigation coefficient profiles. From eqn. (17), it can be seen that the bioirrigation intensity is species-specific, because C_b and C, and therefore the ratio of (C_b - C)/(C₀ - C), depend on the reactivity and diffusion coefficient of the individual solute species.

To assess the significance of concentration changes in burrows, the steady state concentration profile within a burrow can be estimated. Balancing the net delivery of the solute to a given depth by flushing of a vertical burrow with radial exchange across the burrow wall gives

where Q is the volumetric flux of water flushing through the burrow. Q is obtained either by direct measurement or is estimated iteratively using the depth-integrated irrigation coefficient profile. Advective water velocity in the burrow is related to the mass flux of water by u = Q/(πr₁²). If the average concentration of the solute in the bulk sediment remains approximately constant (e.g., in the case of O₂ which is rapidly depleted near the SWI, so that C_avg ≈ 0 over most of the depth range of interest), then the radial derivative in eqn. (18) can be expressed according to eqn. (5), or

Solving eqn. (19) for the boundary condition C_b(0) = C₀ gives

which shows, as expected, that deviation from perfectly flushed burrows is favored by narrow burrows (small r₁), slow water exchange (small u) and high solute reactivity (small L _i ).

The flushing efficiency is related to the residence time of water in the burrows, τ_res = (V_b/Q) where V_b is the burrow volume. If we assume, for simplicity, that α is constant with depth then

Therefore, it is possible to obtain estimates of irrigation coefficients based on flushing frequencies (τ_res^-1), burrow thickness (r₁) and abundance (r₂). Although the specific form of eqn. (21) is only valid in the limited case of α constant with depth, nonetheless, this equation highlights the potential use of data on burrow flushing frequencies of benthic organisms to constrain bioirrigation intensities.

Results and discussion

Model organisms

The burrows of a variety of macrofauna, especially those inhabiting intertidal and shallow subtidal environments, have been characterized by resin casting.[34] Resin casts are obtained by pouring a dense epoxy directly into sediments or mesocosms; after hardening, casts are removed and cleaned. This technique is widely used to study burrow morphology and to obtain estimates of burrow surface areas and volumes. However, it can be difficult to obtain casts of deep-burrowing macrofauna[35] and to distinguish burrows of different types of macrofauna. For example, initial resin cast studies of the echiuran worm Maxmuelleria lankesteri suggested that burrows have a single entrance at the surface. [36, 37] However, subsequent studies of sediment intake and ejection have demonstrated the presence of a second connection to the surface. The second entrance was likely not found in initial resin cast studies because the presence of the large echiuran worms within inhabited burrows blocked the flow of resin [35] and because second openings are often occluded by sediment[38].

X-radiography of sediment slabs from box cores has also been used to estimate burrow shapes and burrow wall surface areas. [32, 39, 40] Burrows of organisms in deeper environments have been studied using underwater television and remote photography, [36, 41, 43] and through mesocosm experiments conducted with organisms collected from sediments[44, 45].

It should be noted that portions of the burrow structure of some organisms may not be irrigated. Therefore, relying solely on resin casting or X-radiography to determine burrow surface areas of exchange may lead to overestimates of irrigation intensity, particularly for deep burrows. This highlights the importance of obtaining detailed information concerning the burrowing and irrigating activies of individual organisms in constructing more representative, ecology-based bioirrigation models.

Decapoda: Thalassinidae: Callianassa subterranea

Thalassinid shrimp are extremely common, inhabiting intertidal, shallow subtidal and possibly deep-sea sedimentary environments. [46, 47] Their burrows, which are used by the organisms for shelter, feeding and reproduction, are deep (sometimes exceeding 2 m depth) and typically have quite complex morphologies. Burrow morphologies may vary for a single species, depending on the sedimentary environment. For example, Callianassa subterranea burrow morphologies are much more complex in sandy [44, 45] than in muddy sediments [46, 48].

Mesocosm studies of C. subterranea collected from sandy North Sea sediments suggest that these shrimp construct burrows with 5.9 ± 1.6 openings at the sediment/water interface connecting to a lattice of nearly horizontal chambers.[45] This complex morphology is approximated in this study by treating each burrow as a combination of three intersecting U-shaped burrows (Table 1). This results in a distinct subsurface maximum in burrow wall surface area at approximately 10 cm depth (Fig. 3A).

Table 1 Parametrization of burrows. Cumulative distribution functions are Gaussian. Ve = vertical burrow parameter; UY = upper branches of Y-shaped burrow parameter; SY = stem of U-shaped burrow parameter; I = 45° inclined burrow parameter; U = U-shaped burrow parameter; Y = Y-shaped burrow parameter, L = L-shaped burrow parameter

Figure 3

Burrow wall surface area per unit sediment volume, S_v, as a function of depth for burrow network simulations of the Thalanassid shrimp Callianassa subterranea using (A) U-shaped burrows (North Sea site parameters), (B) U- and Y-shaped burrows (North Sea & Davies Reef site parameters) or (C) L-shaped and vertical + 45° inclined shapes (Scottish Loch site parameters). Data used for all calculations are shown in Table 1.

However, C. subterranea burrows do not always resemble intersecting U-shapes. Nickell and Atkinson report several C. subterranea burrow casts from muddy Scottish Loch sediments with single, nearly vertical shafts extending downward from a lattice of horizontal tunnels. [48] Tudhope and Scoffin found even more complex C. subterranea burrow networks that typically had 5 or more downward-inclined chambers radiating from a central U-shaped burrow. [49] To account for the presence of the deep shafts found in these studies, C. subterranea burrows are modeled using a combination of 50% U-shaped and 50% Y-shaped burrows (Fig. 3B). The addition of Y-shaped burrows representing deeply penetrating shafts results in the persistence of significant burrow surface area to depths of greater than 80 cm, compared to less than 25 cm in simulations using only U-shaped burrows (Fig. 3A,B). However, the deeper portions of Y-shaped Callianassa burrows may be irrigated only infrequently, so that without consideration of flushing frequency the model may lead to overestimates of the bioirrigation intensity at depth.

In muddy subtidal sediments from the west coast of Scotland, Atkinson and Nash found only one C. subterranea burrow, out of 17 casts, that had more than a single inhalent shaft opening to the surface.[46] However, Nickell and Atkinson found that most burrows from this environment had an additional, much smaller exhalent shaft that opened or partially opened at the sediment surface, perhaps facilitating flow of water through the burrows.[48] Reported organism densities in these sediments are much lower than estimates from North Sea sites, and burrow depths are considerably greater. To examine the potential differences in surface area resulting from the extremes of reported burrow morphologies, these Scottish Loch burrows are modeled here assuming only a single connection to the surface, by using 50% L-shaped and 50% vertical/45° inclined burrow shapes. This results in a distinct increase in burrow wall surface area at approximately 10 cm, which is not dissimilar from the subsurface maxima resulting from simulations using U- and Y-shaped burrows (Fig. 3C). Thus, C. subterranea is likely to promote deep bioirrigation in all of these settings, with a subsurface maximum in burrow surface area at depths of 10 cm or greater.

Echiura: Bonellidae: Maxmuelleria lankesteri

The echiuran worm Maxmuelleria lankesteri has been recognized as an important bioturbator because, unlike many other burrowing macrofauna, it builds extremely deep, temporally persistent burrows. Single burrows have been observed to remain in place for at least 15 months and are suspected to persist in a single location for years. [37, 38] Although M. lankesteri is found from the west coast of Scotland to the Skagerrak, its burrows have been studied particularly in Irish Sea sediments, where the organism is thought to promote deep, rapid mixing of radionuclide contaminants into the sediment [36, 37, 43].

There has been some controversy surrounding the morphology of M. lankesteri burrows. [35, 38] In this study, simulations of M. lankesteri burrows were made using U-shaped burrows, single shaft (L-shaped, inclined L-shaped, 45° inclined, vertical/45° inclined), or mixtures of U-shaped and single shaft burrows. Surface areas as a function of depth are similar for all morphologies; results for mixtures of U-shaped and single shaft burrows are shown in Fig. 4. The simulated burrow surface area profile suggests that M. lankesteri promotes bioirrigation to depths of up to 1 m. However, as seen in Fig. 4, a very high variability in irrigation intensity can be expected below 30 cm.

Figure 4

S_V as a function of depth for burrow network simulations of the Echiuran worm Maxmuelleria lankesteri. Parameters are given in Table 1.

Nereididacea: Nereididae: Nereis diversicolor

The polychaete worm Nereis diversicolor inhabits estuarine, intertidal sediments throughout Europe. [40, 50–52] Very high population densities (2000–4000 individuals m^-2) are common, suggesting that in many estuaries these polychaetes contribute significant bioirrigation intensity. Burrow morphologies have been studied in laboratory mesocosms and field settings using resin casting and X-radiography. [40, 53] Although burrows have frequently been described as having a simple U-shape, Davey has shown that in mesocosms N. diversicolor burrows initiate in a U-shape, but soon exhibit more complex morphologies.[40] After several hours, burrows are extended to form a Y-shape, and after several days have complex morphologies including multiple openings at the sediment/water interface. The morphology of N. diversicolor burrows in natural sediments may depend on animal population densities; U-shaped burrows of N. virens are constructed at low animal densities, but in dense populations burrows may be I-, L- or Y-shaped.[54] Although Gerino and Stora describe N. diversicolor burrows in their study as U-shaped,[53] X-radiographs of these burrows shown in their Fig. 1 suggest that they are Y-shaped. X-radiographs of burrows shown in Davey and Gerino and Stora are used to estimate the depths and horizontal extent of burrows (Table 1). [40, 53] The resulting surface area depth profile agrees closely with burrow wall surface areas measured by Gerino and Stora using resin casts of 4 individual N. diversicolor burrows from a mesocosm study (Fig. 5).

Figure 5

Comparison of S_v for four individual Nereis diversicolor from Gerino and Stora [53] (solid blue line) and from this study (dotted red line). Error bars show one standard deviation calculated from stochastic network simulations.

Decapoda: Ocypodidae: Uca minax, Uca pugnax, and Uca pugilator

Fiddler crabs, including the sand fiddler (Uca pugilator), the mud fiddler (Ucapugnax) and the brackish water fiddler (Uca minax) are among the most numerous and conspicuous of burrowing macrofauna in intertidal mangrove and saltmarsh sediments. [55–58] Estimates of fiddler crab densities in a salt-marsh at Sapelo Island, GA are as high as 205 ± 46 m^-2 and burrows have been reported to penetrate to depths of up to 65 cm. [57, 59] However, unlike other burrowing macrofauna, fiddler crabs do not permanently inhabit or actively irrigate their burrows. In fact, some fiddler crabs plug their burrow entrances to prevent flooding during tidal innundation. [55, 60, 61] Fiddler crab burrows are nonetheless likely to contribute to solute transport, because of the abundance of abandoned burrows that are passively flushed by tidal innundation. However, most studies of fiddler crabs report only aperture or organism densities, making an assessment of flushed burrow density difficult.

Burrow network simulations suggest that U. minax burrows, because of their large size (Table 1), will contribute most to irrigation, especially at depths greater than 30 cm (Fig. 6C). Simulated burrow wall surface area profiles for U. pugilator and U. pugnax both exhibit a subsurface maximum at approximately 10 cm depth (Fig. 6A,B), with few burrows extending deeper than 20–25 cm into the sediment.

Figure 6

S_v as a function of depth for fiddler crabs including (A) Uca pugilator, (B) Uca pugnax, and (C) Uca minax. Parameters are given in Table 1.

Decapoda: Grapsidae: Sesarma reticulatum and Eurytium limosum

The marsh crab Sesarma reticulatum occurs in a variety of intertidal marsh environments.[55–57, 62] Reported population densities range from 30 m^-2 in levee marshes with tall-form Spartina alterniflora to just 1 m^-2 in unvegetated creek banks. [55] Burrows are extensive and, especially near the sediment surface, have complex morphologies. Typical burrows have several entrances attaching a shallowly sloping initial tunnel several cm long to the surface. This initial tunnel network is connected to a single vertical shaft of 2–5 cm diameter that descends into the sediment to depths averaging 13–30 cm, although depths of up to 75 cm have been reported. [56, 62]S. reticulatum burrow networks were simulated in greatly simplified form as single vertical shafts (Table 1). This results in a burrow surface area profile that decays gradually with depth (Fig. 7A).

Figure 7

S_v as a function of depth for the mud crabs (A) Sesarma reticulatum and (B) Eurytium limosum. Parameters are given in Table 1.

Another crab species found frequently in both vegetated and unvegetated intertidal saltmarshes is Eurytium limosum.[55, 57]E. limosum burrows are typically composed of two or more shallow tunnels that extend laterally for 60–70 cm; a single inclined shaft joins the long, shallow tunnels and descends to 20–30 cm depth to create a broadly Y-shaped burrow[57]. E. limosum burrows were simulated in this study by treating each burrow as a combination of two shallow, intersecting L-shaped and one deep 45°-inclined burrows. The presence of just a few E. limosum burrows greatly enhances burrow surface areas in the upper few centimeters of the sediment column (Fig. 7B).

Enteropreusta: Spengelidae: Schizocardium sp

Laboratory mesocosms have been used to assess the role of the funnel-feeding acorn worm Schizocardium sp. in promoting bioirrigation of shallow, near-shore sediments in St. Louis Bay, Mississippi Sound. [32] X-radiographs of mesocosm sediment slabs suggest that Schizocadium sp. build approximately U- or V-shaped burrows penetrating to a maximum of 9 cm depth. Burrow networks were simulated in this study assuming that all burrows are approximately V-shaped. The resulting surface area profiles decrease continuously with depth (Fig. 8A).

Figure 8

(A) S_v as a function of depth for the acorn worm Schizocardium sp. with mean animal density of 500 m^-2. (B) Burrow densities for populations of 100 m^-2 (blue), 311 m^-2 (red), 356 m^-2 (black) and 800 m^-2 (green). Lines are burrow densities from this study; symbols represent number of active burrows reported from X-radiography of mesocosm slabs.[32] Error bars represent ±1 standard deviation in burrow density.

X-radiography of sediment slabs has been used to directly determine Schizocardium burrow densities as a function of depth in four mesocosms. [32] The expected depth-dependence of Schizocardium burrow density, given the animal densities in the various mesocosms (100, 311, 356, 800 m^-2), is calculated with the stochastic model with parameters given in Table 1 (Fig. 8B). Because total animal densities are known and therefore are kept invariant in the stochastic model, standard deviations shown in Fig. 8B are relatively small. In all four cases, burrow density profiles calculated using the stochastic model are in very good agreement with reported densities from X-radiography of sediment slabs.

Dry Tortugas, Florida, USA

Macrofaunal consortia

The consortium of bioturbating macrofauna present in shallow-water carbonate reef sediments at Dry Tortugas, FL has been described by D'Andrea and Lopez.[63] Of the deeply (>4 cm) bioturbating organisms, D'Andrea and Lopez suggest that the polychaete worm Notomastus sp. (density = 113.2 m^-2) and the burrowing shrimp Callianassa sp. (density = 40.8 m ^-2) are dominant.[63]Notomastus, a deep deposit feeding capitellid polychaete, was commonly found at depths greater than 15 cm, most often occurring at 25–30 cm depth[63].

To simulate burrow networks at Dry Tortugas, only the dominant deep bioturbators Notomastus and Callianassa are considered here. Notomastus burrows are assumed to be similar in morphology to those of the polychaete Nereis diversicolor, with penetration depths of 20 ± 2.5 cm for Notomastus. In fact, this is a great simplification; Notomastus burrows have been reported to resemble complex spirals (D'Andrea and Lopez, personal communication), thus, Notomastus burrows may contribute a greater surface area of exchange than is considered here. Callianassa sp. are assumed to build burrows with morphology similar to those reported for Great Barrier Reef Sediments[49] and are therefore modeled using the U- and Y-shaped Callianassa burrow parameters given in Table 1. The presence of these two organisms alone leads to a burrow surface area profile with three distinct regions; in the upper 6 cm surface areas are high and essentially constant with depth, from 6–12 cm they decrease rapidly and then more gradually in the depth interval 12 to 90 cm (Fig. 9A).

Figure 9

(A) S_v as a function of depth in shallow carbonate sediment at Dry Tortugas, FL, USA, from simulations using a consortium of the polychaete worm Notomastus sp. and the shrimp Callianassa sp. (B) Bioirrigation coefficients (s^-1) as a function of depth for Dry Tortugas, FL derived by fitting an early diagenetic model to measured pore water profiles (Furukawa et al.;[32] red solid lines), and estimated from the burrow network in this model. Blue line indicates α_O2 calculated using L = 2.6 mm from Furukawa et al.,[15] with eqn. (5), assuming C_b = C₀. Green line indicates corrected for C_b using eqn. (17) and (20). (C) Bioirrigation coefficients (s^-1) as a function of depth from Furukawa et al.[32] (red lines) and estimated from the burrow network in this model. Blue line indicates α_SO4-2 calculated using sulfate concentration and reduction rate profiles from Furukawa et al.;[32] green line indicates α_SO4-2 corrected for C_b, using eqn. (17) and (20). Orange solid line indicates bioirrigation coefficient profile calculated using eqn. (16), assuming zero order kinetics for sulfate reduction. Dotted lines indicate irrigation coefficient profiles calculated using ±1 standard deviation of the stochastic model burrow densities and surface areas.

Bioirrigation coefficients

From the surface area profiles shown in Fig. 9A, bioirrigation coefficients are calculated using eqn. (6) and (10). D _i is calculated using porosity and molecular diffusion coefficient data given in Furukawa et al. [15] and r₁ is calculated using eqn. (7) with surface areas and burrow densities obtained from the stochastic network simulations. The O₂ penetration depth at the SWI (L_O2) is set to 2.6 mm, based on vertical microelectrode profiles measured by Furukawa et al.,[15] and is used to calculate according to eqn. (10); values vary from 4.8 to 6.4 mm. At the SWI, this results in predicted irrigation coefficients that are approximately twice as high as the values estimated by Furukawa et al. from diagenetic modeling of chemical data[15] (Fig. 9B), perhaps reflecting the use of L_O2 measured at the SWI, rather than at the BSI. The predicted irrigation coefficients also exhibit a less pronounced decrease with depth than the coefficients obtained by Furukawa et al.[32] The higher irrigation coefficients estimated in this study may be due, at least in part, to imperfect burrow flushing. Eqn. (17) and (20) are used to assess this possibility. This correction leads to lower irrigation coefficients, with mean values within a factor of 2 of estimates from Furukawa et al.[15] Nonethless, the results suggest that O₂ gradients measured at the SWI may not yield accurate estimates of radial diffusion length scales at the BSI at depths greater than a few mm.

Bioirrigation coefficients are also calculated by combining the stochastic burrow network simulations with sulfate concentration and reduction rate profiles taken from Furukawa et al.[15] The resulting irrigation coefficients predicted using eqn. (6) and (12) are only slightly larger than those predicted by Furukawa et al.[15] and exhibit a similar depth-dependence (Fig. 9C). Because sulfate is a less reactive solute than O₂, correction of this profile for inefficient flushing using eqn. (17) and (20) has negligible effect (Fig. 9C). Both the corrected and uncorrected sulfate bioirrigation coefficient values are significantly smaller than those calculated for O₂. This is primarily because of the much larger radial diffusion length scales calculated for SO₄^-2. Sulfate reactive length scales are also large relative to those obtained from O₂ microprofiles, because they are calculated from bulk sulfate reduction rates that decrease rapidly with depth, reflecting changing reactivity with depth.

Bioirrigation coefficient profiles are also calculated here completely independently from data in Furukawa et al.,[15] by assuming zero order kinetics for sulfate reduction and using eqn. (16). In this case, the irrigation coefficients, which depend only on the burrow sizes and spacings, are approximately a factor of 3 lower than those calculated by Furukawa et al.[15] right at the SWI (Fig. 9C). Below a few mm, however, good agreement between the two independent approaches is observed.

Sapelo island, GA, USA: unvegetated creek bank

Macrofaunal consortia

Intertidal saltmarshes at Sapelo Island, GA are populated by diverse, abundant burrowing macrofauna.[55, 57, 59, 64–66] Basan and Frey report 1040 burrow apertures m^-2 in an unvegetated creek bank at Sapelo Island,[57] and suggest that the majority of these apertures are due to polychaete worms including Nereis succinea and Heteromastus filiformus. Other burrows are attributed to fiddler crabs, especially Uca pugilator, burrowing shrimp including Upogebia affinis and predatory mud crabs such as Panopeus herbsti, which share the burrows of Sesarma reticulatum.[56] Teal also inventoried burrowing macrofauna at an unvegetated creek bank on Sapelo Island, GA and found Uca pugilator (52 m^-2), Uca pugnax (13 m^-2), Sesarma reticulatum (1 m^-2) and Eurytium limsosum (6 m^-2)[55].

A simulated burrow network for an unvegetated creek bank at Sapelo Island, GA was completed using U. pugilator, U. pugnax and E. limosum densities taken directly from Teal,[55] while for S. reticulatum and P. herbsti averages of data given for S. reticulatum by Teal and for P. herbsti by Basan and Frey are used.[55, 57]Upogebia affinis shrimp burrows are assumed to be similar in morphology to Callianassa shrimp. A polychaete density of 450 ± 50 m^-2 is derived by assuming a total aperture density of ~1000 m^-2 of which ~100 m^-2 are due to organisms other than polychaetes.

Burrow wall surface areas for polychaetes, shrimp and fiddler crabs only are highest near the sediment-water interface and decay gradually to a depth of approximately 30 cm (Fig. 10A). The addition of E. limosum/P. herbsti burrows results in much higher burrow wall surface areas in the upper few cm of sediment and a more rapid decay of the burrow surface areas with depth (Fig. 10B).

Figure 10

(A) S_v as a function of depth in an unvegetated saltmarsh creekbank at Sapelo Island, GA, USA including polychaete worms, fiddler crabs and the mudcrabs Sesarma reticulatum and Panopeus herbsterii. (B) S_v (cm²/m³) with all of the organisms included in (A) as well as Eurytium limosum. (C) Bioirrigation coefficients (s^-1) as a function of depth calculated using inverse modeling of chemical data (red solid line)[68] and estimated from the burrow network in this model using L = 1.0 mm with eqn. (6), assuming C_b = C₀ (blue solid line) and with correction for depletion of O₂ within burrows (green line). (D) α_SO4-2 bioirrigation coefficients calculated using the stochastic network model with measured sulfate reduction rates, without correction for C_b (blue line) and with correction for C_b (green line). Orange line indicates bioirrigation coefficient profile calculated with the stochastic network, assuming zero order kinetics for sulfate reduction (eqn. 16). Dotted lines indicate irrigation coefficients calculated using ±1 standard deviation of the burrow density and surface area profiles.

Bioirrigation coefficients

As for the Dry Tortugas site, the stochastic model-derived surface areas shown in Fig. 10B are used with eqn. (6) and (10) to derive a bioirrigation coefficient profile. A porosity of 76%[67] is used to obtain , and r₁ is calculated using eqn. (7). Using a measured O₂ flux of 68–89 mmol O₂ m^-2 d^-1 yields O₂ penetration depth () values between approximately 0.2 mm and 1.0 mm.[72] Given that the marsh sediments are periodically exposed to the atmosphere, an value of 1.0 mm is used in this study, which results in values that vary from 4.1 to 17.0 mm. The resulting bioirrigation coefficient profile, shown in Fig. 10C, is compared to that of Meile et al.[68] who used an inverse modeling approach to obtain a bioirrigation coefficient profiles from sulfate concentration and reduction rate measurements. In the upper 5 cm of sediment, bioirrigation coefficients from this study are approximately a factor of 2 greater than those of Meile et al.,[68] and the bioirrigation coefficient profile from this study does not decay as rapidly with depth. Eqn. (17) and (20) are used to correct the profile for the effects of imperfect burrow flushing. The resulting profile still suggests more intense irrigation throughout the sediment than the profile of Meile et al. [68].

The lower irrigation values given by Meile et al.[68] may reflect differences in solute-specific irrigation coefficients. Therefore, irrigation coefficients for sulfate are calculated here using sulfate concentration and reduction rate profiles [67, 68] and eqn. (6) and (12). The resulting irrigation coefficient profile (Fig. 10D) is in excellent agreement with the profile given by Meile et al. [68] Correction of the stochastically-derived bioirrigation coefficients for the effects of imperfect flushing using eqn. (17) and (20) does not significantly change the profile.

Bioirrigation coefficients are also calculated assuming zero order kinetics for sulfate reduction. The resulting profile, which depends only on the burrow sizes and spacings, and which is completely independent from the sulfate concentration and reduction rate profiles used by Meile et al., [68] is also shown in Fig. 10D. The resulting bioirrigation coefficients are somewhat smaller than those of Meile et al.[68] Overall, the results for the creek bank sites indicate that it is possible to obtain meaningful irrigation coefficients from the stochastically simulated burrow network.

Sapelo island, GA, USA: vegetated ponded mash

Macrofaunal consortia

The abundance and composition of macrofaunal communities present in saltmarsh environments shows significant spatial zonation, with community composition depending on length of tidal innundation, presence or absence of vegetation and sediment composition.[57] Thus, areas of the marsh with dense, diverse macrofaunal populations, such as the unvegetated creek bank discussed above, are likely to have deeper, more intense bioirrigation than regions of the marsh with fewer organisms. For example, ponded marsh regions are much more sparsely populated with macrofauna than creek bank sites. At Sapelo Island, GA, Basan and Frey report ~415 burrow apertures m^-2 in the ponded marsh,[57] compared to 1040 m^-2 at the unvegetated creek bank. Teal found only two types of crab in the ponded marsh, the fiddler crab Uca pugnax (27 ± 7 m^-2) and the mud crab Eurytium limosum (4 ± 1 m^-2).[55] All other burrow apertures (~280 m^-2) are likely due to the polychaete worm Nereis succinea. Basan and Frey report that burrow depths are typically shorter in the densely vegetated ponded marsh,[57] and that burrow shapes are typically less complex. In simulations, Uca pugnax burrows at the ponded marsh were assumed to reach a mean depth of only 15 cm, compared to 20 cm at the creek bank. The lower density of organisms inhabiting the ponded marsh leads to lower surface areas as a function of depth compared to the creek bank site (Fig. 11A).

Figure 11

(A) S_v as a function of depth for the ponded marsh (solid line) sites at Sapelo Island, GA with dotted lines indicating ±1 standard deviation of the surface area. (B) Bioirrigation coefficients (s^-1) as a function of depth at Sapelo Island calculated for dissolved O₂ at the creek bank (green) or ponded marsh (black) with correction for solute depletion within burrows and without correction for solute depletion at the creek bank (blue) and ponded marsh (purple). (C) Bioirrigation coefficients (s^-1) as a function of depth at Sapelo Island calculated for dissolved SO₄^-2 with correction for solute depletion within burrows at the creek bank (green) or ponded marsh (black) and without correction for solute depletion at the creek bank (blue) and ponded marsh (purple).

Bioirrigation coefficients

Bioirrigation profiles at the ponded marsh are calculated for both dissolved O₂ and SO₄^-2. O₂ irrigation coefficient profiles are calculated using a vertical diffusion length scale of 1.0 mm, as for the creek bank site. Diffusion coefficients were corrected using the porosity of 85% reported by Kostka et al. [67] Dissolved SO₄^-2 irrigation coefficient profiles are calculated using measured sulfate concentration and reduction rate profiles.[67, 69] Correction for inefficient burrow flushing results in slightly lower O₂ irrigation coefficients, but has no significant impact on SO₄^-2 irrigation coefficients (Fig. 11B,C). The intensity of irrigation, especially as indicated by the O₂ irrigation coefficient profiles is somewhat less than at the creek bank site. The depth of irrigation is also somewhat shallower at the ponded marsh, because of the lower density of deep-burrowing infauna and because infauna in ponded marsh areas tend to build shallower burrows than in creek bank sediments.[57] Because of shallower irrigation, the ponded marsh sediments have a more compressed redox stratification, with more reducing conditions closer to the sediment surface, compared to the more intensely irrigated creek bank sediments[69].

Comparison of diagenetic and stochastic irrigation models

Model results obtained at Dry Tortugas, FL and Sapelo Island, GA demonstrate that the ecologically-based stochastic model yields bioirrigation coefficients that are comparable to those obtained independently via inverse modeling or by use of multicomponent early diagenetic reactive transport models. The latter do not require information concerning the benthic organisms responsible for irrigation, rather, bioirrigation coefficients are inferred from model simulations fitting measured chemical and/or reaction rate profiles (Meile et al., 2001). In contrast, the stochastic model presented here explicitly uses ecological data to calculate irrigation intensities. Pore water profiles are not per se required, although the radial diffusive length scale (L) across the burrow-sediment interface must be specified. As few direct measurements of L are currently available, L values are constrained here using other available chemical data.

An important strength of the stochastic model is that uncertainties associated with calculated bioirrigation intensities can be evaluated. The uncertainties derive from the probability density functions describing the burrow network, which in turn are based on statistical information about the ecological variables (animal density, burrow size and geometry). Thus, the model-calculated uncertainties on the bioirrigation coefficients reflect the natural variation in abundance and burrowing activity of the benthic macrofauna.

The stochastic model is not intended to replace early reactive transport models of early diagenesis. Rather, it is an additional tool to constrain irrigation coefficients, which can then be used in early diagenetic modeling. Independent estimates of transport parameters significantly enhance the reliability of reactive transport calculations applied to the complex biogeochemical dynamics of sediments.

Conclusion

Aquatic sediments are typically redox stratified, with increasingly reducing conditions encountered at depth. However, burrow networks built by macroinfauna create a complex 3-dimensional patchwork of relatively more oxidized or reduced sediment compartments. The resulting increase in interfacial area separating zones with disparate redox characteristics enhances diffusive exchange and creates ecological and geochemical microzones that affect sediment biogeochemical dynamics.

In this study, the spatial distribution of distinct sediment zones due to burrowing macroinfauna are quantified by simulating 3D burrow networks for a variety of model organisms and for consortia of organisms in shallow carbonate and intertidal saltmarsh sediments. The conceptual approach provides a link between the ecological characteristics of burrowing macroinfauna and the resulting biogeochemical consequence on solute transport and biogeochemical cycling. The impact on diffusional exchange of solutes between burrow water and bulk sediment is addressed, and quantitative estimates for nonlocal transport parameters are proposed. This is of particular importance, because a quantitative description of biogeochemical dynamics in sediments, for example based on the interpretation of chemical profiles, requires the quantification of transport intensities.

The approach outlined here allows an independent, ecologically-based, assessment of the depth-dependence of bioirrigation coefficient profiles. Burrow network simulations of even relatively simple model organisms are found to result in depth-dependent burrow wall surface areas that are neither constant nor exponentially decreasing over a given interval, as is commonly assumed in bioirrigation models. Network simulations of the burrowing shrimp Callianassa subterranea, the echiuran worm Maxmuelleria lankesteri, and the fiddler crab Uca pugilator all result in burrow wall surface areas characterized by a subsurface maximum. These maxima are similar to subsurface maxima in bioirrigation intensity reported in chemically-based studies that have not imposed a priori restrictions on the depth-dependence of bioirrigation intensity.[27, 68] Furthermore, model simulations demonstrate that relatively rare organisms may have a disproportionate influence on solute transport in sediments. For example, burrow wall surface area profiles calculated for an unvegetated saltmarsh creek bank site at Sapelo Island, GA with and without the mud crab E. limosum are quite distinct (Fig. 10A,B), in spite of the fact that E. limosum, is responsible for only 6 out of ~1000 apertures m^-2.

Burrow network simulations provide estimates of burrow surface areas as a function of depth, which are used to derive nonlocal bioirrigation coefficients. In carbonate reef sediments at Dry Tortugas, O₂ based bioirrigation coefficients calculated using the stochastic model are somewhat larger than values found through diagenetic modeling.[15] However, if the irrigation coefficients are corrected for inefficient flushing, using an estimate of the time-averaged value of C_b, the agreement is better.

Sulfate bioirrigation coefficient profiles calculated using the stochastic model for an unvegetated intertidal saltmarsh creek bank site are in good agreement with nonlocal bioirrigation coefficient profiles determined by Meile et al.[68] using an inverse modeling approach. Calculated irrigation coefficients for an adjacent ponded marsh site are somewhat lower. This lower mixing intensity may contribute to the more compressed vertical redox zonation observed in chemical pore water profiles at the ponded marsh compared to the creek bank[69].

To go beyond the model presented here, further information is required regarding the depth-dependence of solute concentrations within burrows, as well as the temporal variation of solute concentrations within burrows. Nonetheless, the current model provides an extremely useful link between benthic macrofaunal ecology and a quantitative description of chemical transport by bioirrigation.

Note

†Presented during the ACS Division of Geochemistry symposium 'Biogeochemical Consequences of Dynamic Interactions Between Benthic Fauna, Microbes and Aquatic Sediments', San Diego, April 2001.