Burrowing in marine muds by crack propagation: kinematics and forces

Burrowing animals dominate marine sediments that constitute 70% of Earth's surface. Burrowers are considered ecosystem engineers, significantly altering their environments and microbial activities (cf. Meysman et al., 2006). In spite of their numerical and biomass dominance, however, they have been underrepresented in the literature on animal locomotion. Progress has been limited in two ways: observationally (marine sediments being literally `clear as mud') and conceptually, by a lack of understanding of sediment mechanics on spatial and temporal scales relevant to burrowing macrofauna.

Newton's third law is commonly applied in biological fluid dynamics, for example to estimate drag forces on a body by measuring the opposite force(exerted by the animal on the fluid), i.e. the rate of momentum extraction from the fluid. The Navier–Stokes equations permit this estimate through explicit fluid parameters of density and dynamic viscosity. In this subfield of continuum mechanics it would be inconceivable to study animal locomotion without considering these material properties. That no study has estimated forces and work of burrowing by reference to explicit parameters of the solid elastic (or viscoelastic) continuum through which burrowers move and against which metabolically fueled forces operate signals a rudimentary state of understanding.

Forces exerted by both marine and terrestrial burrowers have been measured in various ways, although no methods, to our knowledge, have yet explicitly considered sediment mechanics. Coelomic (internal) pressure has been measured with a cannula through the body wall (e.g. Ansell and Trueman, 1968; Seymour, 1969; Seymour, 1971; Trueman and Foster-Smith,1976; Hunter and Elder,1989). In many of these studies, burrowers were close to rigid walls so that behaviors could be observed (e.g. Seymour, 1971; Trueman and Foster-Smith,1976; Hunter and Elder,1989). Close proximity to a rigid interface increases stiffness of the sediment against which the animal exerts forces(Dorgan et al., 2006). Measured internal pressures from which forces are calculated are higher near walls than animals would exert in the natural environment to achieve the same burrow/crack opening (body thickness). In addition to cannulae, Hunter and Elder (Hunter and Elder, 1989)used an isometric force transducer attached by hook and thread to the tail of the worm to measure `tail-pulling' force as an empirical measure of the work needed to overcome friction. Trevor(Trevor, 1978) used a diaphragm to separate sediment from water connected to a pressure transducer to measure forces exerted against the diaphragm by the anteriors of burrowing worms. Force transducers have also been used to measure anterior and radial forces (Quillin, 2000). Diaphragms and force transducers, however, offer resistances different from natural soils or sediments. Similarly, internal pressures of worms crawling on the surface or moving in water (e.g. Seymour, 1969) are not representative of internal pressures of worms burrowing in natural sediments(Dorgan et al., 2006).

Bubbles in muddy sediments create disk-shaped cracks that grow and permit bubble rise by fracture. Their growth and aspect ratios have been modeled using linear elastic fracture mechanics (LEFM) theory, indicating that muddy sediments behave in a linear elastic manner on these small time and space scales (and under force magnitudes) typical of burrowing behavior(Johnson et al., 2002; Boudreau et al., 2005). LEFM theory considers three material properties for two-dimensional (2-D) problems in which the material is isotropic and the crack is loaded only in mode I(opening or uniaxial tension): elastic modulus (E), critical stress intensity factor or fracture toughness (K_Ic), and Poisson's ratio (ν). E is a measure of stiffness and the constant that relates stress, σ (force/area), and strain, ϵ(elongation/original length), as σ=Eϵ in a linear elastic,isotropic material undergoing uniaxial deformation. Poisson's ratio(dimensionless) is the negative of the constant of proportionality between longitudinal and transverse strain under uniaxial stress. An incompressible material such as gelatin has a Poisson's ratio of 0.5, whereas lower values of Poisson's ratio indicate higher compressibility. Mode I stress intensity, K_I, is the coefficient of the dominant term in the series expansion of the stress field at a crack tip under mode I loading, and is used to compare stresses at cracks of varying configurations under varying loading conditions. The crack propagates when K_I exceeds the critical value, K_Ic, the fracture toughness. Another way of stating the fracture criterion is that when energy release rate(G) exceeds resistance of the material (R), the crack grows. If energy release rate increases as the crack grows, growth is unstable (e.g. shattered glass), whereas if energy release rate decreases as the crack grows,growth is stable and stops when energy release rate falls below material resistance. Relevant fracture mechanics have been recently reviewed (cf. Dorgan et al., 2006; Mach et al., 2007) and can be found in textbooks (e.g. Anderson,1995).

Adhesive and cohesive forces of the mucopolymer matrix holding grains together (resulting in elastic behavior) dominate mechanics of muds, whereas clean, monodisperse sands are granular materials for which the weights of individual grains are more important in determining contacts(Dorgan et al., 2006). Burrowing mechanics differ between the two media, and only muds are considered here.

The burrow around the polychaete Nereis virens Sars is a planar crack (shaped like a tongue depressor) and extends laterally away from the worm, with elastic rebound of the medium compressing the worm dorsoventrally(Dorgan et al., 2005). The most recently produced crack segment is an oblate hemispheroid that propagates when the worm everts its pharynx to exert dorsoventral forces against crack walls. These forces concentrate stress at the crack tip, producing fracture when stress intensity exceeds the critical value, K_Ic (cf. Anderson, 1995).

Using photoelastic stress analysis (cf. Harris, 1978; Full et al., 1995; Dorgan et al., 2005), a technique once widely employed by engineers to observe stress patterns (cf. Durelli and Riley, 1965), we measured forces exerted by Nereis virens burrowing in gelatin, a clear analog for muddy sediments. We quantified differences in relevant material properties between gelatin and muddy sediment and used finite element modeling to calculate forces exerted by the polychaete N. virens in muddy sediments from forces measured in gelatin. Finite element analysis has largely replaced photoelastic stress analysis to evaluate stress distributions in engineering applications yet, as shown here, the photoelastic method remains useful in an experimental context. This paper gives the first estimates of forces exerted during burrowing in natural sediments that explicitly consider the mechanical properties of sediments.

To enhance visibility of stress fields, large (>5 g wet mass) specimens of Nereis virens Sars were selected. Animals were obtained from Harbor Bait (Edgecomb, ME, USA) and kept in containers of mud under flowing seawater until use.

Critical stress intensity factor (K_Ic) and elastic modulus (E) for seawater–gelatin and mud have been measured(Johnson et al., 2002). For mud from Cole Harbour, NS, USA, K_Ic is 280–490 Pa m^1/2 and for gelatin, K_Ic is 50–220 Pa m^1/2. They measured E of sediment as approximately 1.4×10⁵ Pa and of gelatin as 1.5–10×10³ Pa(Johnson et al., 2002; Boudreau et al., 2005). Considerable variation exists, and sediment has higher E and K_Ic than gelatin, but the ratios K_Ic/E coincide approximately based on ranges of values for the parameters, and bubbles observed in both media have predictable and similar aspect ratios (Boudreau et al.,2005). The relevance of sediment mechanics to burrowers has been recently reviewed (Dorgan et al.,2006).

Because aspect ratios of bubbles depend on the K_Ic/E ratio of the medium, it seems reasonable to assume similar dependence for aspect ratios and extensions of animals'crack-shaped burrows. One difference between gelatin and muddy sediment is a greater loss of stored elastic potential energy in sediment than in gelatin,resulting in a lower relative elastic restoring force and less dorsoventral compression of the body in mud (see Results). However, the total forces exerted by animals in sediment are higher because stiffness is higher, so elastic restoring forces may be comparable in the two media. We have observed difficulty by Sternaspis burrowing in gelatin (appearing compressed and often failing to burrow; K.M.D. and P.A.J., unpublished observations), but Nereis seems to burrow without obvious difficulty. We are working to develop a better analog, but have not yet formulated a better, isosmotic,non-toxic material with the transparency necessary for photoelastic stress analysis.

For gelatin, E was calculated using Eqn 2 from displacements resulting from test tubes of known mass resting upright on the surface of aquaria used for experiments (aquaria being too large for the Vitrodyne tester). We also tested smaller containers of gelatin with the Vitrodyne tester (as described for sediment) for comparison.

Gelatin, a birefringent material, was placed between two circularly polarizing filters of opposing polarizations, with a light source on one end and a camera on the other (Harris,1978; Full et al.,1995; Dorgan et al.,2005). When filters initially are lined up, no light passes. Stress in the material re-orients light in directions of maximal and minimal stress; that re-oriented light passes through the filter, showing up as a lighted region. Area and intensity of the resulting light and dark patterns are related to the state of stress in the material(Harris, 1978).

The experimental setup (Fig. 1) comprised a Just Normlicht (Weilheim/Teck, Germany) Smartlight 5000 photographic light table covered first with a green color filter (Rosco CalClor #4430, Rosco Laboratories Inc., Stamford, CT, USA), and then a right-handed, circular polarizing filter (3M HNCP 37% R.H. S-10×0.030 in from Edmond Optics, Barrington, NJ, USA). In front of the filter was a 20.8-l glass aquarium of double-strength (2×) seawater–gelatin (28.35 g l^–1). A CCD videocamera (Basler A622f, Exton, PA, USA)(camera 1) with 6× close-focus zoom lens (Edmund Optics #52-274,Barrington, NJ, USA) was sited opposite the light table to record in lateral view images of the worm (defined as the y–z plane) at 7.5 frames s^–1. On the lens were a 52-mm, green (061) color filter, then a 52 mm, left-handed (standard), circular polarizing filter. Because the Rosco color filter is a gelatin sheet, it was placed between the light and the polarizer; placing it between the polarizers would show stress in the color filter, interfering with images. An identical camera (camera 2)and lens at 90° recorded the dorsal (or ventral) view of the worm (defined as the x–z plane) at 3.75 frames s^–1. A smaller light table (Porta-trace/Gagne 10×12 in) was sited opposite camera 2 and was covered by a magenta filter (Rosco CalColor #4790), then a right-handed, circular polarizing filter. On the camera lens were a 52 mm,magenta (CC30M) filter, then a left-handed (standard), circular polarizing filter. The color filters partitioned the light spectrum between the two cameras, avoiding interference from the orthogonal light source. Cameras were attached to separate computers, and videos were recorded digitally and analyzed using LabView software (version 7.1.1, National Instruments, Austin,TX, USA).

We used circularly rather than linearly polarizing filters (cf. Full et al., 1995) because circular polarizers show a larger lighted region (cf. Sharples, 1981). Circular polarizers combine a linear polarizer with a quarter-wave retarder(Fig. 1). Our images show a single patch of light resulting from a force, and the area of the light field is proportional to the force.

Experiments were conducted in a cold room at 11°C. The only light sources were the photographic light tables, which were completely covered with the filters. Food-grade gelatin(www.bulkfoods.com)was made with seawater boiled to reduce viable bacteria and set overnight(12–36 h) in the cold room.

Each tank of gelatin was calibrated using test tubes resting upright on the surface with volumes of water added (0.2–4.0 ml) as known mass. Aperture and zoom of the camera lens were set before each calibration and kept constant throughout experiments. Images were captured for each mass, and force regressed against number of pixels lighter than a threshold value for each tank (Fig. 2; see Appendix). Thresholds were chosen for each tank to be as high as possible while lighting an area of ⩾300 pixels on each side of the worm's everted pharynx. Photoelastic fringes, lines separating light from dark regions, indicate contours in the stress field. Because we worked with small stresses, we saw only primary fringes. There is a linear relationship between force and area of the primary compression fringe (Harris,1978; Full et al.,1995), but in our system that linear relationship holds for only a limited range of areas, with smaller stresses showing quadratic size changes. The observed area of the primary compression fringe is a 2-D projection of a 3-D stress field, and the order of the relationship is expected from mechanics to fall between linear (resulting from a large, distributed load) and quadratic (resulting from a point load). The range of pixel areas around worms falls mostly below the linear range, and a second-order polynomial was fitted(see Appendix). We first digitally removed the light intensity gradient near the surface of the gelatin due to stress caused by evaporative shrinkage by subtracting an image taken after the test tube was removed, then adding a background pixel value. The background pixel value was an average of pixel values in an undisturbed dark region in the middle of the aquarium and was constant for all images used in a calibration. Images were analyzed using LabView (National Instruments).

Test tube width (0.01135 m) was chosen to make the diameter of contact of the test tube with the gelatin surface (0.0053–0.0079 m, depending on mass) close to the mean width and length of the worms' everted pharynges(0.0067 and 0.0068 m, respectively). The curved bottom of the test tube gave a distinct patch of light more similar in shape to the patches around the worms than did flat, rigid cylinders in preliminary tests. Test-tube width provided a scale in the images to attach a length to the pixel dimension. Experimental validation of the calibration method is presented in the Appendix.

After calibration, a crack was initiated perpendicular to the view of camera 1 with forceps, and a worm placed within. If worms pulled back out of the crack, they were gently replaced until they started burrowing. The macro lenses used have fixed focal lengths, so cameras were moved to keep worm distance from the camera constant.

We used only videos of the lateral view of the worm in which the plane of the crack was perpendicular to the lens of camera 1 (parallel to the x–z plane cf. Fig. 1). Body twisting, diagnosed as a widening of the crack tip anterior to the worm, caused video rejection. If the worm moved toward or away from camera 1 (greater than 30° from vertical as observed with camera 2),the light regions were too large (because some body stress lined up with the stress from the everted pharynx), and videos were rejected. To restrict wall effects, we used only videos in which worms were >0.05 m away from all aquarium walls. At 0.05 m, the effect of the wall on K_I is approximately 10%, calculated by comparison to the exact solution for an edge crack in a finite-width plate.

Segments of video (3.75 frames s^–1) fitting the above criteria and the additional criterion that the camera was not moved (worm moved parallel to both camera planes) were used to measure frequency of pharynx eversion, distance moved between pharynx eversions, and resultant velocity. LabView was used to measure coordinates of the anterior tip of the worm, the crack tip, and width of the pharynx.

Frame grabs around each pharynx eversion fitting the above criteria were converted to binary images, and pixels above threshold counted in LabView. Compressional stresses on dorsal and ventral sides of the pharynx were analyzed separately because a separate force was applied to each crack wall. Thresholded images showed, on each side, primary compression fringes produced by the pharynx, tension at the crack tip and, in some images, internal body pressure just posterior of the everted pharynx(Fig. 3). The `kidney-bean'shaped tensile stress fringe agrees qualitatively with stress patterns at crack tips under stresses perpendicular to the plane of the crack. To measure force exerted by the everted pharynx, we included only pixels in the region of compressive stress. In some cases, compressive and tensile stress fringes were indistinguishable, and these images were rejected. Forces were calculated from number of pixels through the calibration curve.

The projected planar area of the pharynx could only be measured directly from the corresponding images from camera 2 when the worms were moving straight down in the aquarium (in an x–z plane). For consistency, we instead measured pharynx width from camera 2 and length from camera 1 and approximated planar area as an ellipse. Internal pressure of the pharynx was calculated as the average of the dorsal and ventral forces divided by the planar area of the everted pharynx.

Internal body pressure was calculated from lighted areas of high stress around dilated segments of the body (that could be seen moving anteriorly as a peristaltic wave). Light regions were smaller and less intense than around the pharynx, so lower thresholds were used on the same images from the calibration. The area exerting the force was calculated as the planar area of one peristaltic wavelength, the product of the width of the worm (from camera 2) and the length of the region exerting the stress, measured as the distance along the worm between two regions (measured from the middle of the contact area of each). Internal pressure was calculated as the force divided by this area. Most of the visible stress was in the anterior of the body; very little stress was visible in the posterior. Stress was rarely visible in the absence of a clear peristaltic wave. Because we measured only visible stress, our measurements are closer to stress maxima than to averages.

Internal body pressure and pharynx pressure were used as inputs for finite element modeling of worm shape in gelatin with the program franc2d (Cornell Fracture Group, Cornell, NY, USA). Franc2d is a two-dimensional, finite element modeling program designed for fracture that calculates both displacements resulting from applied stresses and stress intensity factors at crack tips. Because a worm burrows in 3-D and the model is only 2-D, we developed two models for two different views of the worm: lateral(y–z plane) and anterior (x–yplane). Both models ran in plane strain mode, which assumes that the thickness of the material was large (i.e. not a thin plate) and that all loads,geometric parameters and solution fields were independent of the through-thickness coordinate. They can only approximate actual 3-D configuration of the burrowing worm, but results discussed later show that a 2-D analysis provides important insights into burrowing mechanics.

We first developed the lateral-view model, a rectangular geometry with the dimensions of the front (camera 1) view of the aquarium. Preliminary models of only the worm's pharynx greatly underestimated displacements, but modeling the worm's body (approximately the length of worms in gelatin during experiments)as well as the pharynx produced displacements much closer to those observed. An edge crack starting at the top surface extending half-way down the aquarium(0.1 m) represents the worm's burrow. Positions of the bottom and sides of the rectangle were fixed, assuming that the gelatin in the aquarium is stuck to the glass walls and does not move, as observed in experiments. Average stress exerted by the pharynx was applied to the crack walls from the crack tip to a point 0.00725 m behind the crack tip (the length of the crack wall needed to get a final displaced pharynx length of 0.00667 m; see Results). Calculated internal body pressure was applied along the rest of the crack wall to 0.014 m from the top surface (model 1). Preliminary model trials showed that applying body stress all the way to the top lifted the top surface of the modeled tank much higher than observed in the aquarium of gelatin, likely because gravity is not included in the model. Because observed body stresses were visible only in the worms' anterior regions, the model was also run with linearly decreasing body stress (model 2), from the measured maximum body stress to the maximum stress that the worms could exert without producing visible stresses(intercept of the calibration curve). Because of the plane strain assumption,the model ignores lateral crack edges, potentially overestimating displacements.

To evaluate the importance of lateral crack edges as constraints on displacements, we developed an anterior model. A rectangular geometry was used with dimensions of the top of the aquarium. An interior crack of length 0.0066 m (the width of the worm) represents the anterior view of the worm's burrow(model thickness equals length of worm). Bottom, top and sides of the rectangle were fixed, as they represent the four sides of the aquarium. Calculated internal body pressure was applied to the crack walls. With stress kept constant over the original crack length, crack tips were extended byΔ a and resultant maximum displacements (δ_max)calculated.

Stress intensity factors were calculated for the lateral model. The median value of three different methods of calculating K_I using franc2d was compared to the critical stress intensity factor for gelatin,50–220 Pa m^1/2 (Johnson et al., 2002).

Stiffness in gelatin (E_gel) was then increased to that of sediment (E_sed), and Poisson's ratio (ν) was decreased from 0.45 to 0.3. Although gelatin is incompressible and has νclose to 0.5, we used a slightly lower value (0.45) to accommodate displacements in a plane strain model. Poisson's ratio for soils varies widely(from 0 to 0.5) depending on the soil type, confining pressure and saturation state (Lade, 2001). Saturated soils are incompressible, but observations of burrowing animals indicate that on these small spatial scales, forces result in compression of the solid-phase sediment by dewatering of the polymer–sediment matrix (explaining the presence of permanent burrows in sediments). Although Poisson's ratio for incompressible, saturated sediments is technically 0.5, the linear elastic model does not take into account dewatering on small spatial scales and longer time scales that result in small, permanent deformations (see Results). Using a lower value of Poisson's ratio is not technically correct, but is a reasonable way to approximate nonlinearities using a simple linear model and worked well in our model validation experiments (see Appendix). We also calculated forces for a Poisson's ratio of 0.45 for comparison. We multiplied stresses measured in gelatin by E_sed/E_gel to calculate approximate stresses that the worms need to apply in natural sediments to have the same body shape. Increased stiffness requires proportionally higher stresses to obtain similar displacements. These stresses were input into the model and resulted in larger displacements because Poisson's ratio had been decreased(increasing compressibility). We then reduced the stresses until displacements matched observed body thicknesses in gelatin. Pharyngeal stresses were then converted back to forces by multiplying by planar pharynx area. Stress intensity factors were calculated and compared to critical stress intensity factors for sediment, 280–490 Pa m^1/2(Johnson et al., 2002). This modeling method was tested for a known system, and the results are presented in the Appendix.

Elastic moduli of muddy sediments (calculated from Eqn 2) were variable: E_sed=(27±10)×10³ Pa (mean ±s.d.; N=8). Although we used large cores (0.15 m diameter), there may have been some wall effects, with overestimation of E. Removal of the core, however, caused sediment to settle considerably and expand radially. Because sediment was disturbed by the removal of the core and absence of constraint from surrounding sediments is more artificial than the constraining walls of the core, we did not use data from unconstrained cores. Values of E_sed decreased when the core was removed (from 1.57 to 1.05×10⁴ Pa in one sample), as expected. We also removed surface layers from the unconstrained core and found that E increased with depth (from 1.05 to 1.61×10⁴ Pa at 0.01 m depth and 2.04×10⁴ Pa at 0.03 m in one sample), but stayed near the range of variability observed across cores.

Elastic modulus for gelatin was E_gel=(1.9±0.3)×10³ Pa (mean± s.d.; N=4). We also measured E in smaller (0.10 m diameter, 0.10 m deep) containers of gelatin with the Vitrodyne V-1000 tester and found higher values of E (=5.6×10³ Pa) that we attribute to wall effects from the necessarily smaller container.

Representative loading–unloading curves(Fig. 4) for mud and gelatin show elastic behavior for both materials. Gelatin is clearly linearly elastic,whereas the sediment curve shows lower resilience as well as a small plastic deformation following initial loading. Subsequent loadings remain elastic,however, and show slightly less loss of stored energy.

Mean distance traveled between pharynx eversions was 0.0073±0.0018 m(±s.d.; N=6). Average time between pharynx eversions for worms that moved without stopping was 9.5±3.8 s (N=6), for an average velocity of (8.7±3.5)×10^–4 m s^–1 (N=6) (Table 1). Worms sometimes stopped moving, but would often continue again after the tail was gently touched. Most worms exhibited consistent frequencies of pharynx eversions and distances traveled, and this behavior was similar for most worms observed (including the many samples eliminated because forces could not be measured). However, we observed one worm that traveled at approximately twice the velocity of average worms and did so with very few pharynx eversions. Unfortunately, this worm was not oriented with the plane of the crack in line with camera 1, so forces and body width could not be measured.

A typical worm moves forward, extending the crack, then begins pharynx eversion while continuing forward and extending the crack. The crack tip does not extend beyond the anterior of the pharynx as it is being everted and moving in the anterior direction (Fig. 5). The pharynx reaches its most anterior point before full eversion (Fig. 5A). As the pharynx moves back and reaches its maximum width, the crack tip is visible. Between eversions, the worm moves its head from side to side within the plane of the crack, extending it laterally with the palps (see Movie 1 in supplementary material). Antennae often extend to and probe the crack tip,tracing its edge.

Forces exerted by the everted pharynx were measured for six worms, five of which exerted forces between 0.014 and 0.016 N(Table 2). The sixth worm exerted much smaller forces, 0.007 N, but pharynx width, length and thickness were smaller than in other worms. Because we were looking for a maximum force exerted and it appeared that the sixth worm was not fully everting the pharynx(jaws not visible), we rejected that data point. Mean force exerted by the other five worms was 0.015±0.001 N. Average pharynx area was(3.7±1.1)×10^–5 m²(Table 2). The six worms used had mean mass of 8.2±3.9 g.

Displacements in the lateral model (Fig. 6) are much closer to the shape (thickness) of the worm than displacements in the anterior-view model(Fig. 7A), which are over an order of magnitude smaller than the thickness of the worm. Extending crack tips by Δa while leaving applied stresses constant increases modeled thickness, but maximum displacement reaches an asymptote after increasing by a factor of only 2–2.5, still much smaller than observed displacements (Fig. 7B).

The two dimensions included in the 2-D lateral-view model are length and dorsoventral thickness. Changes in the lateral direction, specifically the lateral crack edges, are assumed to be unimportant in constraining the shape of the worm. We tested that assumption using the anterior model, in which lateral edges are present and maximum displacement increased as the lateral constraint was removed by extending crack width (cf. Fig. 7B). The observed distance between the lateral edge of the worm's body and the lateral edge of the crack(visible in the right panels of Fig. 5B) was 0.0047 m. At Δa=0.0047 m, maximum displacement was 0.70 times that as Δa approached infinity, the assumption in the lateral model (cf. Fig. 7B). Because the anterior model suggests that lateral edges were constraining thickness of the worm's body, we ran lateral models with constant(model 1) and linearly decreasing (model 2) body stresses reduced by this constraint factor of 0.70 (models 3 and 4). Decreasing body stress by the constraint factor resulted in decreased body displacements and slightly decreased pharynx displacements (Fig. 8, Table 3).

Thicknesses of everted pharynges for observed worms are larger than the modeled displacements in the lateral-view model. We ran another model,extending the length over which the pharynx stress was applied from the average to the maximum observed pharynx length (model 5). Modeled pharynx thickness was larger than in the other models, but still smaller than observed thicknesses (Fig. 8).

Stress intensity factors (K_I) for the five models(Table 3) ranged from 57 to 64 Pa m^1/2, within the range of critical stress intensity factors for gelatin, 50–220 Pa m^1/2(Johnson et al., 2002). This result supports the use of the simplified lateral model and, more importantly,validates the mechanism of burrowing by crack propagation as mechanically feasible. Worms are capable of generating stress intensities that can drive crack growth.

Modest discrepancy between observed and (lateral) modeled displacements is reasonable considering the simplified, 2-D, plane–strain approximation of the 3-D observations and considerable variability in material properties included in the model. In the model, stresses are applied perpendicular to the crack, whereas stresses applied by the worm as the pharynx everts likely have radial components. Stresses in other directions could not be included in the model, which is run in a single step. The worm pushes the pharynx forward,which could explain why the observed pharynges are thicker than model predictions. Also, in the model, constant stress is applied over the entire pharynx area, but the force may be applied over a smaller area than the measured planar pharynx area and is unlikely to be evenly distributed over that smaller area. More focused modeled stresses would result in a more pronounced displacement peak. Variation in measured E for gelatin could explain some of the discrepancy, and there are likely other sources of error in approximating a 3-D process with a 2-D model. One specific example is the use of Poisson's ratio (ν) of 0.45 instead of 0.5, the value for an incompressible material. An incompressible material in plane strain with fixed walls cannot deform because there is nowhere for the deformed material to go;setting ν=0.5 unsurprisingly caused the program to crash. Using a lower value of Poisson's ratio suggests, incorrectly, that the gelatin is compressible, but is necessary for the 2-D model to work at all; this compressibility in the model roughly equates with gel expansion upward in the real aquaria.

An alternative explanation is that error lies not in the modeling assumptions, but in measured forces input into the model. It is possible that the force measured from compression fringes may be slightly underestimated because the compressive stress fields might be affected by the tensile stress field at the crack tip. We believe that this error, if present, is very small because the compressive stress regions are small and usually appear to be clearly isolated from the tensile stress fields. It seems much more likely that the discrepancy comes from the simplifying assumptions of the model rather than experimental error.

Given the potential sources of error in the model, displacements in all five models are reasonably close to observed values, and modeled stress intensity factors fall within the range of K_Ic for gelatin, adding support to model representations of experimental results. Unfortunately, the results from the five models(Fig. 8, Table 3) were not close enough to observed worm shapes, which themselves showed extensive variability, to determine which model is most appropriate. All five use the same value of pharynx stress; although body stress does affect pharynx thickness,differences are not critical in using the model to convert forces exerted in gelatin to those in natural sediments. If we assume that the discrepancy between modeled and observed worm shapes results from simplifying assumptions in the model, it is reasonable to further assume that the same discrepancy will occur in modeling a worm in natural sediments. It follows that any of the models could be used to convert forces from gelatin to natural sediments as long as the shape of the worm in the sediment model matches that in the gelatin model.

We selected model 2, with unconstrained, linearly decreasing body stress,to calculate force exerted by the pharynx in natural sediments. We increased E from E_gel (=1.9 kPa) to E_sed (=27 kPa), decreased ν from 0.45 to 0.3, then calculated new stresses by multiplying stresses in the gelatin model by E_sed/E_gel. Pharyngeal stress was increased from 378 to 5305 Pa and body stress from 92 (linearly decreasing to 60) to 1291 Pa (linearly decreasing to 842). Applying the new stresses to the model resulted in larger displacements than in the gelatin model because of increased compressibility (from the change in Poisson's ratio), and stresses were scaled down to obtain a similar worm shape(Fig. 9). Stresses applied in the final model were 4408 Pa for the pharynx and 1073 Pa (linearly decreasing to 700) for the body. We then multiplied pharyngeal stress by average pharynx area, 3.7×10^–5 m², to calculate force exerted to propagate a crack in natural sediments, 0.16 N.

We followed similar methods to obtain the force for ±1 s.d. of the elastic modulus for sediments, (2.7±1.0)×10⁴ Pa. The forces needed to match the modeled displacements for E=(1.7 and 3.7)×10⁴ Pa are 0.10 and 0.22 N, respectively. We also calculated the force required for E=1.39×10⁵ Pa to be 0.83 N, as measured in deeper sediments(Johnson et al., 2002). We followed the same procedure without changing Poisson's ratio, and calculated a force of 0.20 N for E=2.7×10⁴ Pa, within the range of variability of E.

Burrowing forces and mechanics depend on mechanical properties of sediments, not only on each of E and K_Ic, but also their ratio. The force in sediment depends directly on E,assuming that the shape of the worm remains constant. An increase in E requires a larger force to obtain the same displacement. However,exerting a larger force makes the crack propagate more easily by increasing the stress intensity factor, K_I, above the critical value, K_Ic. Because we assume constant displacement and stress exerted depends on E, resulting stress intensity factors depend on the value of E. For E=(1.7, 2.7 and 3.7)×10⁴ Pa, K_I=(4.7, 7.4 and 10.2)×10² Pa m^1/2, respectively, compared to K_Ic of 2.8–4.9×10² Pa m^1/2 for sediments (Johnson et al., 2002). Stress intensity factors within or above the range of critical values for sediments support the mechanism of burrowing by crack propagation in sediments as well as in gelatin.

The larger values of K_I are intriguing, although further research is needed to determine if they are significantly higher than K_Ic. One possible explanation is that our assumption of constant body shape is inaccurate; if worms are flatter in muddy sediments,less stress would be exerted, resulting in a lower K_I. More likely, our assumption of linear behavior of sediments is inaccurate over longer intervals (relevant to body stresses). Applied body stress may be lost to frictional dissipation, creep or plastic deformation rather than amplified at the crack tip.

Infrasound measurements have shown that activities of individual animals in muddy sediments can be detected as pressure waves up to 0.5 m away(Wethey and Woodin, 2005). Such pressure signals, which also depend on the mechanical properties of sediments, could be used as an independent measure of burrowing forces if mechanics and heterogeneity of sediments were well enough known.

Burrowing forces further depend on body size: radial forces exerted by earthworms scale with (body mass)^0.43(Quillin, 2000). We did not consider allometry of burrowing forces but expect smaller worms to exert smaller forces than the large nereidids used in this study.

Crack growth can be characterized as stable or unstable; stable crack growth is associated with displacement-driven fracture (e.g. a wedge driven into a piece of wood, creating a crack opening as thick as the wedge), whereas unstable crack growth is associated with load-driven fracture (e.g. the two sides of the splintered wood are pulled with constant force, and as soon as the force exceeds a critical value, the wood breaks in two)(Anderson, 1995). Crack extension as the worm moves forward agrees with descriptions in the fracture mechanics literature of stable, wedge-driven crack propagation. The entire body acts as a long wedge, and pharyngeal eversion thickens the wedge at the tip of the crack where it has greatest effect. Poor fit of the anterior model to the observed thickness of the worm in gelatin, as well as poor fit of a preliminary lateral model in which only the pharynx was modeled, suggest that the body acts as a (relatively) unconstrained wedge. The worm moves from side to side in the burrow between pharyngeal eversions, extending the crack laterally. This behavior is likely important in removing this lateral constraint, allowing greater displacements of the medium along the body with less stress than if crack walls were closer to the worm. Worms are able to move forward without eversion, both over short distances in periods between pharyngeal eversions and over longer distances. This mode is consistent with observations of other worms that burrow by crack propagation without pharyngeal eversion [e.g. Leitoscoloplos fragilis (Polychaeta;Orbiniidae), Heteromastus sp. (Polychaeta; Capitellidae), and Saccoglossus kowalewskii (Hemichordata); K.M.D. and P.A.J.,unpublished observations]. Stable, wedge-driven fracture by worms is also consistent with the wedge shape of hard-bodied burrowers such as clams and urchins.

Energetic costs of burrowing have been considered much higher than for other forms of locomotion (Trevor,1978; Hunter and Elder,1989). Previous estimates, however, were also calculated from force measurements rather than measured directly through calorimetry or indirectly through oxygen consumption. External energy was calculated from measured forces, distances over which those forces were applied(force×distance=work), and the animal's velocity. Multiplying external energy by an energetic efficiency yields a net cost of transport. In previous studies, forces were often measured against rigid walls, overestimating them[e.g. 0.56 and 0.67 N radial forces exerted by Polyphysia crassa and Priapulus caudatus, respectively(Hunter and Elder, 1989),compared to our calculation for Nereis of 0.16 N]. In addition,because mechanics of burrowing and of the medium were not understood, the distances over which the forces were applied were assumed to be distances the animal moved (Trevor, 1978)rather than much smaller distances perpendicular to the direction of motion.

We originally thought that the primary problem with the use of external energy to calculate net cost of transport for burrowing, once forces and distances were accurately measured, was that we did not know energetic efficiency. However, our modeling work suggests that measuring external energy use is much more complicated. Force exerted by the everted pharynx and distance perpendicular to direction of movement can be multiplied to estimate work, but modeling shows that neglecting body stress underestimates displacements and results in stress intensity factors below critical for fracture. Furthermore, the worm must also exert propulsive forces parallel to the direction of motion in order to move forward in the burrow. Although it is possible to calculate total external work from applied stresses using franc2d(by integrating along the length of the crack), the model assumes a linear elastic material without loss of stored energy or creep over time. Sediment does behave linearly on short time scales, but the body of the worm applies pressure longer than does the pharynx, and the linear assumption is unlikely to hold over those longer periods. Attempts to link forces exerted by burrowing animals to energetic costs of burrowing (e.g. Hunter and Elder, 1989; Trevor, 1978) are complicated by the nonlinear and poorly understood potential energy storage and loss behaviors of sediments that need additional modeling and measurement attention at spatial and temporal scales of burrowing.

List of symbols and abbreviations

 
• a

  half-length of crack (m)

 
• a_c

  half-length of bubble (i.e. crack) at fracture (m)

 
• B

  body stress

 
• C

  compressive stress

 
• d

  diameter of cylinder (m)

 
• E

  elastic modulus/Young's modulus (Pa)

 
• E_gel

  elastic modulus for gelatin (Pa)

 
• E_sed

  elastic modulus for sediment (Pa)

 
• F

  force (N)

 
• G

  energy release rate (potential energy available for fracture) (J m^–2)

 
• I_p

  influence factor for shape of a rigid plate

 
• K_I

  stress intensity factor (Pa m^1/2)

 
• K_Ic

  critical stress intensity factor (Pa m^1/2)

 
• LEFM

  linear elastic fracture mechanics

 
• P

  pressure

 
• R

  resistance to fracture (J m^–2)

 
• T

  tensile stress

 
• Δ

  displacement of cylinder (m)

 
• δc

  half-thickness of bubble (i.e. crack) at fracture (m)

 
• δmax

  maximum displacement (i.e. thickness of worm) (m)

 
• Δa

  change in crack length, crack growth (m)

 
• ϵ

  strain (dimensionless)

 
• ν

  Poisson's ratio (dimensionless)

 
• σ

  stress (Pa)

Because worms burrow in 3-D, we conducted tests to evaluate whether the linear relationship between stress and the area of the primary compression fringe observed in shallow, flat gelatin plates(Harris, 1978; Full et al., 1995) holds for deep, 3-D aquaria and whether a surface calibration is appropriate to measure forces around a worm burrowing within the gelatin. We first calculated regressions between pixel areas and stresses and forces applied to the surface by objects of varied sizes and shapes, to determine both whether the linear relationship between stress and area of pixels held for 3-D stress fields and also whether the geometry of a test tube was appropriate to mimic the everted pharynx. These objects included balloons with varying volumes of water,Play-doh™ (Hasbro, Pawtucket, RI, USA) objects of varying weights and shapes, and test tubes and flat-bottomed cylinders of varying sizes. Using a larger test tube resulted in larger pixel areas for a given stress(Fig. A1A), but plotting force instead of stress gave similar results across test tubes sizes(Fig. A1B). Curved pieces of Play-doh™ and balloons showed similar stress fields to those of test tubes, but flat, rigid objects had smaller patches of stress along the edges with little stress in the middle of the object, very different than the stress fields around worms. Rigid, flat objects create constant displacements in gelatin, exerting high stress at the ends and low stress in the middle (a parabolic stress field), most noticeable by photoelastic stress analysis for light objects exerting small stresses.

We then used finite element modeling to apply simulated stresses over different areas. A 2-D axisymmetric model of half of the gelatin tank (the x–z plane) with stress applied to the surface from the top corner of rotation (corresponding to the center of the aquarium) to a distance simulating the radius of a test tube was used to evaluate the effects of radius and magnitude of stress on areas of primary compression fringes(high stress). Again, the relationship between force and pixel area was much more linear than between stress and pixel area(Fig. A1C,D). Relationships for both calibration tests are linear rather than quadratic because the forces are larger (compare to x values in Fig. 2).

Because worms move within gelatin, exerting forces on two walls, and our calibration involves only one surface, we also compared the test-tube calibration with data from a balloon inflated in the gelatin. The balloon was stretched over a curved wire to form a flat disk and was attached to a syringe and a 1 PSI gauge pressure transducer (Honeywell, Columbus, OH, USA) connected through a data-acquisition device (National Instruments USB 6008) to a computer running LabView (Fig. A2). The balloon was filled with enough water to give a pressure reading close to zero, inserted into the sediment (or gelatin), inflated to known volume, and pressure measured. The balloon was deflated, removed from the gelatin, then re-inflated in air to the same volume to determine the pressure needed to inflate the balloon without resistance. Subtracting this pressure removed the stress needed to stretch the balloon, leaving only the stress applied to the gelatin. This approach assumes that balloon and gelatin stiffnesses are additive, which is approximately valid in this experiment since balloon and gelatin displacements must remain compatible. Forces exerted were calculated by multiplying pressure by the planar area of the balloon in contact with the gelatin. The planar area of the balloon overestimated the area over which stress was applied for small volumes of water that did not fully inflate the balloon. We were unable to quantify the error in the planar area measurements, but this error seems to explain the discrepancy between the regressions using the test tube and balloon(Fig. A3A). Overestimating the area over which stress is applied overestimates the calculated force(=stress×area) for small volumes of water; shifting those points to the left would bring the balloon regression closer to the test-tube regression,although the variation around the balloon regression is still very high. Although the experimental data were less conclusive than we had hoped, they did not reveal any obvious problems with the surface calibration method.

Finite element modeling to evaluate the difference in geometry between the crack tip and the surface of the gelatin was more convincing. Stress was applied in a 2-D plane strain model (again, a simplification of the 3-D state)to the top surface of modeled gelatin and along the wall of the crack near the tip (Fig. A3B). Applying stress of 500 Pa along 0.0089 m starting at the crack tip resulted in a slightly smaller region above the threshold stress than the same stress applied along the same length at the crack surface (42 and 64 pixels,respectively). However, moving the stress 0.002 m back from the crack tip (the distance between the fully everted pharynx and the crack tip; see Results)yielded more similar stress distribution to the surface (58 and 64 pixels,respectively). Reduction in compressive stress fringe size at the crack tip likely resulted from interference by the tensile stress at the crack tip (cf. Fig. 3), visible as an asymmetry in the tensile stress field in Fig. A3Bii compared to Fig. A3Biii. This interference would be overestimated in the model compared to the experiment because observed stress fringes are narrower and extend farther from the surface than modeled stresses. This difference suggests that stresses in the experiment are less uniform than modeled stresses because the curved shape of the test tube(and worm pharynx) exerts more stress in the center than at the edges and/or the model ignores effects of displacements on shape of the stress field.

We performed the calibration with several different threshold values to calculate forces exerted. Calculated force increased with threshold value at low thresholds, reaching an asymptote at higher threshold values (data not shown). This increase is likely due to interference between compressive stress around the pharynx and tensile stress at the crack tip. Using a higher threshold value yields a smaller area of the primary compression fringe that is better resolved from the tensile stress field.

The slope of the quadratic calibration curve was higher for smaller numbers of pixels; using too high a threshold increased error in converting from pixels to force. However, using low thresholds underestimated force because of influence of the tensile stress field. To restrict error, we used the highest threshold (close to the asymptote) that resulted in areas of at least 300 pixels.

To test for appropriate displacements under applied stresses, we modeled a control system with known stress and displacement. A balloon attached to a pressure transducer and inflated with a syringe(Fig. A2) was inserted through a vertical slot in the side of a 16.5×21 cm container of sediment that had been fully mixed and allowed to settle for 2 weeks. The same methods as for the calibration validation were used, and the experiment was repeated in both sediment and gelatin with two balloon sizes. Volumes of the combined balloon, wire and added water were measured.

Balloons were modeled in franc2d as an internal crack in an axisymmetric model of the container of sediment (or gelatin). Measured pressures were applied to the crack and volume of displaced sediment was compared to measured volume of the balloon and water. In both media, model results approximated actual volumes (Table A1).

The authors would like to thank Eric Martin and Bob Lad from the University of Maine Laboratory for Surface Science and Technology (LASST) for use of the Vitrodyne tester, Dennis Hill of Harbor Bait, Edgecomb, ME, USA for providing worms, and Bernie Boudreau, Bruce Johnson, Eric Landis and Larry Mayer for helpful discussions and comments on the manuscript. This research was funded by an NSF Graduate Research Fellowship and an NDSEG fellowship to K.M.D. and by ONR Grant No. N00014-02-1-0658 to L. Mayer.