The beach system is located at the northern side of the back barrier island of Spiekeroog. The nearshore and offshore zone of Spiekeroog is exposed to the hydrodynamics of the open sea and strongly affected by tides, winds and wave driven processes (Badewien et al., 2009; Lettmann et al., 2009; Grashorn et al., 2015; Holinde et al., 2015; Schulz et al., 2015). Recent modeling studies suggest that density gradients are among the major drivers of long-term coastal dynamics (Burchard et al., 2008; Burchard and Badewien, 2015). With an average tidal range of 2.7 m at its western end and 2.9 m in the front of the mainland, Spiekeroog represents a standard mesotidal system (Hild et al., 1999). Characteristic are shoreface connected ridges in front of the island coastline (Markert et al., 2015). Those ridges running from southwest to northeast in water depths from 4 to 8 m and from southeast to northwest in water depths of 9 to 20 m. Inflowing high tide waters (flood), as well as outflowing low tide waters (ebb) are guessed to pass these reefs alongshore and offshore. The effects of current and wave driven processes result in different seasonal influences on the system (Dobrynin et al., 2010; Grashorn et al., 2015). During sampling time in May 2014, the beach system was in a seasonal transitional stage and can, thus, be classified as intermediate stage (McLachlan et al., 1996; Short, 1996; Defeo and McLachlan, 2005). The study site at the near and offshore zone of Spiekeroog island included 17 sampling stations, which are located along a 3.4 km long transect line perpendicular to the shore line (Figure 1). The distance between the sampling stations was 0.2 km. The topography of the sampled transect shows a smooth increase in depth with no remarkable structures along the transect line. The beach slope is 0.3°.

All measurements were conducted at calm weather conditions during the 21st to the 23rd May 2014 and obtained with the RB “Otzum” and the RV “SENCKENBERG.” Continuous hydroacoustic and current measurements were conducted on 13 north-south and 7 west-east orientated transects. These transects included the transect, where samples for macrofauna and sediment analyses were taken (Station 1–17). At the starting and end point of each transects, a profiling CTD was used for determining water temperature, salinity and depth.

At each of the 17 stations, sediment samples were sampled from the 0.1 m² van Veen grabs taken for macrofauna identification. After a photographic documentation of the sediment surface and a macroscopic description, subsamples of the sediment surface (the upper 2 cm) were taken from the grabs for grain size analyses.

In the laboratory, the sediment samples were desalted by means of fresh water drainage over 24 h at semipermeable hoses. To separate the mud fraction (< 0.063 mm), the desalted sample was separated using a 0.063 mm mesh screen. The coarse fraction was then dried and dry-sieved over a 2 mm mesh screen to separate gravel and shell debris (coarser 2 mm). To determine the content of shell debris in the sand fraction (>0.063 to < 2 mm), differences of the weight content were registered before and after the elimination of carbonats using hydrochloric acid. The weight content of the shell debris of both grain size fractions (gravel and sand) were simply summed up to get the total amount of shell debris content of the entire samples. The sand fraction was analyzed in a settling tube by means of hydraulic grain size calculation on the base settling velocities measurements by a Macro-Granometer^TM (Flemming and Thum, 1978; Brezina, 1979). Following the Wentworth (1922) grain size scale the higher resolved hydraulic grainsize distribution were normalized to 0.25 phi intervals in order to the mechanical sieved coarse fractions.

A parametric echosounder coupled with a dual-frequency SideScanSonar system (SES2000 standard plus INNOMAR^TMRostock) was used to survey the 13 north-south orientated and 7 west-east orientated transects covering the transect sampled for macrofauna (Figure 2). It was surveyed by means of the parametric echosounder and dual-frequency SideScanSonar system SES2000 standard plus (INNOMAR^TM Rostock). SES2000 sub bottom profiler was running at a primary frequency of ~100 kHz, secondary frequency was 8 to 10 kHz. The SideScanSonar system recorded data on the base of 250 kHz. A swath range of 75 m was selected. Chirp was modulated by linear frequency mode (LFM).

Bathymetry was extracted of the bottom track of the generated SES vertical profiles. High resolution positioning was done with Trimble852 Real Time Kinematic system, with an accuracy of < 5 cm in all 3 directions. The data was further processed and georeferenced by ArcGIS10.22013 Esri.

In a first step, a hierarchical cluster analysis followed by a similarity profile test (SIMPROF; Clarke and Gorley, 2006) was conducted based on square root transformed species abundance data using the statistical PRIMER package, v6 (Clarke and Gorley, 2006) to identify macrofauna communities. Similarities were calculated using the Bray-Curtis coefficient. ANOSIM randomization test (Clarke and Gorley, 2006) was used to test significant spatial differences in the community structure. Discriminating species, that significantly (p < 0.05) separated the different communities, were identified using the similarity percentage routine (SIMPER). The evenness (J') and the Shannon Wiener Index (H') were calculated using the DIVERSE tool.

Also spatial patterns for hydrodynamic, sediment and hydroacoustic parameters of the beach system were analyzed separately by a hierarchical cluster analysis followed by a similarity profile test (SIMPROF; Clarke and Gorley, 2006) based on Euclidean-distance-similarity of normalized environmental data sets of the sediment and the hydrodynamic parameters. For the spatial seafloor characterization, classified acoustic backscatter data and ground truth information of visual and grab data were compiled.

In a second step, the relationships between the macrofauna communities and all environmental parameters combined were analyzed by a joint hierarchical cluster analysis followed by a similarity profile test (SIMPROF; Clarke and Gorley, 2006) based on Euclidean-distance-similarity and performed on the environmental groups. Distance from the shore line was included as categorized variable. Data were further processed and visualized by ArcGIS10.22013 Esri.

Based on this analysis, a spatial zonation of the beach system could be defined. This previous zonation was improved by a principal component ordination analysis (PCO) performed on Euclidean-distance-similarity of environmental parameters using the PERMANOVA+ add on PRIMER v6 (Anderson et al., 2008). The most important environmental parameters were chosen according to Pearson's correlation on the PCO.

A total of 85 taxa were identified. Thirty six taxa were polychaetes, 29 crustaceans, 14 bivalves, 3 echinoderms, and 3 represented other taxonomic groups. Polychaetes were by far the most abundant taxonomic group constituting 90.3% of the total abundance of recorded species, followed by bivalves (5.2%), echinoderms (2.8%) and crustaceans (1.7%).

The multivariate cluster analysis with square root transformed abundance data, based on SIMPROF (p < 0.05), revealed three clusters significantly separated (ANOSIM 0.85; p < 0.0001) and two stations, which were not allocated to a cluster (Figure 3, Table 1). The two non-allocated stations were stations 1 and 2 on the transect line located near to the shoreline. The three clusters and two stations represented five statistically significant different macrofauna communities according to different abundances of the characterizing species Lanice conchilega, Magelona johnstoni, and Spiophanes bombyx.

Communities 1 to 3 (Table 1) revealed mean taxa numbers of 19–21/0.1 m², mean abundances of 38–154/0.1 m², Shannon-Wiener diversities of 1.61–2.11 and evenness of 0.53–0.71. The characterizing species Lanice conchilega, Magelona johnstoni and Spiophanes bombyx occurred in different abundances in the communities, while Bathyporeia elegans was the dominant species in Community 2.

Community 4 (Table 1) revealed higher mean taxa number and mean abundance than Communities 1–3, but Shannon-Wiener diversity and evenness were a bit lower. Lanice conchilega was the dominant species, followed by Spiophanes bombyx and Magelona johnstoni.

Community 5 (Table 1) was characterized by a species-rich benthic fauna and highest mean abundances, while lowest Shannon-Wiener evenness were found. Abundances of Lanice conchilega reached 1702 ind./ 0.1 m².

The normalized hydrodynamic data obtained from the sampling stations for depths of the bottom boundary layer (BBL) and the upper bottom boundary layer (UBBL) were clustered hierarchically based on similarity. Looking at current velocities within the two layers, differences in depth (Δz) and in current velocity (Δu), as well as differences in temperature (Δt), salinity (Δs), and direction of outflowing low tide waters, we identified five hydrodynamic groups (A-E), which were significantly separated by SIMPROF (p < 0.05). The underlying environmental data of the five groups is shown in Table 2.

Current velocity in the UBBL increased continuously along the transect from group A to E from 0.25 to 0.40 m/s, while in the BBL it was highest in group D resulting in low difference in velocity of both layers (Δu), similar to group A (Table 2). Temperature was 1.73°C lower in group E than in group A. Salinity increased from group A to E. In groups A and B, low tide waters flowed in northwest direction, in groups C to E in southwestern direction.

The hierarchical clustering based on similarity of the normalized environmental data of depth, bulk contents of sand, mud and shell debris, and mean sand grain sizes revealed five sediment types (a-e). Within this types only type a was significantly separated by SIMPROF (p < 0.05). The underlying environmental data of the five sediment types is shown in Table 2.

Sandy sediments dominated. Sand contents decreased generally along the transect from 96 (type a) to 90% (type e), while lowest content was found in type c with 84%. Content of shell debris ranged between 1.78 (type a) and 3.53% (type d) (Table 2). Contents of medium and fine sand were highest in type a, but content of medium sand decreased significantly from type b onwards, while content of fine sand increased (Table 2).

Based on backscatter data of the SideScanSonar and the ground truth data of the grab samples, the study site could be subdivided into three different areas along the transect (Table 2).

Area 1 was found within water depths of ~2 to ~3 m at low tide level. It included sampling stations 2 to 4, which backscatter signals were classified as areas of coarser sand with virtually no shell fragments.

Area 2 was found within water depths of ~5 to ~6 m at low tide level. It included sampling stations 7 to 11, which were classified as mixed areas, since there were some muddy areas, as well as fine sand areas with shell fragments and smaller worm tubes of the species Lanice conchilega. Some fecal pellets covered smaller shell fragments and the sand material.

Area 3 was found within water depths of ~8 to ~11 m at low tide level. It included sampling stations 13 to 17. Photographs showed an increase in shells of the bivalve Ensis directus and dense aggregations of Lanice conchilega tubes.

The cluster analysis of all environmental parameters and the macrofauna communities subdivided the sandy beach system of the near- and offshore zone of Spiekeroog island into three zones (Figure 2, Table 3): The first zone included stations 1 to 6 in the inner nearshore zone, the second zone included stations 7 to 12 in the outer nearshore zone, and the third habitat included stations 13 to 17 in the offshore zone.

The principal component ordination analysis (PCO) of the environmental parameters confirmed the separation of the stations within three zones depending on their distance from the shoreline (Figure 3). The outcome states that 61.1% of the spatial structure shown by PCO can be explained by the environmental parameters. The first environmental axis was explained by the hydroacoustic seafloor classification, and by the hydrodynamic parameters (temperature, salinity, Δu, trajectories of outflowing low tide waters). The second environmental axis was associated with the sediment parameters (shell debris content, mud content, fine sand, medium sand).

Within the inner nearshore zone, we found three different sediment types: sediments, which contained high contents of sand; sediments, which contained high contents of fine sand; and sediments with increased mud contents (< 3% to ~ 6%). Several studies already showed the relationship between hydrodynamic mediated sediment composition and spatial distribution of macrofauna communities (e.g., Rachor et al., 2007; Kröncke et al., 2011; Schückel et al., 2013; Armonies et al., 2014; Markert et al., 2015). The physical environment can significantly influence food availability and effect sediment stabilization (Brown and McLachlan, 1990; Short, 1996; Defeo and McLachlan, 2005; Defeo et al., 2009). Small changes in any environmental parameter can form small microhabitats for macrofauna organisms (McLachlan et al., 1993; Defeo and McLachlan, 2005). Since in the highly dynamic environment of sandy beaches habitats are under permanent pressure and can change within short-time scale, macrofauna has to adapt quickly. Interesting in our study was for example the decrease in abundances and the change in species composition in Community B and C, where the sediments were characterized by an increase in fine sand contents probably due to decreasing hydrodynamic energy.

We found three different macrofauna communities in the inner nearshore zone, representing a comparably small distance (~1 km). Communities were characterized by higher diversity indices, while taxa numbers and abundances were lower indicating an intermediate level of disturbance (Clarke and Gorley, 2006). McLachlan et al. (1993) postulates by the swash exclusion hypothesis (SEH) that swash climate, mainly influenced by wave height and beach slope, is the main driver controlling the intertidal sandy beach macrofauna. Defeo et al. (2001) formulated the Habitat Harshness Hypothesis (HHH), which states that in harsh reflective beaches, organisms need to divert more energy toward maintenance, and, therefore, have lower fecundity and higher mortality. However, the decrease in abundances found in the Community B might have resulted from a punctual increase in the “harshness of the environment” reflected also by greater depth. The occurrence of characterizing species such as Bathyporeia elegans and Haustorius arenarius in Community B, which are common in medium sands and high exposition rates to hydrodynamics such as the inner nearshore zone of the nearby island of Norderney (Dörjes, 1976), support a strong influence of hydrodynamics in the exposed inner nearshore areas of sandy beach systems. Higher mud contents and lower depth in Community C revealed the high spatial variability and reflect higher food supply for species such as Lagis koreni, Ensis directus, and Macoma balthica probably caused by nutrient discharge mediated growth of phytoplankton and microphytobenthos (Beck et al., 2017).

In the outer nearshore zone the environmental parameters generally revealed more stable conditions ongoing with changes in the community composition and distribution. We found a slight increase in the current velocities of the BBL and the UBBL. The trajectories of the current velocities changed and implied alongshore and offshore transports. Based on the assumption that the transition point, characterized by the change in the direction of outflowing low tide water, marks the end of the breaker and surf zone (Short, 2012; Carter, 2013), influences of wave action on the sea bottom in the outer nearshore zone should be minimal. This fits to our findings of the sediment distribution in the outer nearshore zone. Sediments were quite homogenous and mainly consisted of fine sand material with low mud contents. With regard to the adapted Autecological Hypotheses of McLachlan (1990), the abundances and taxa numbers of macrofauna communities should now increase due to the decreased “harshness of the environment” in the outer nearshore zone; biological interactions should become more important (Defeo and McLachlan, 2005; McLachlan and Dorvlo, 2005). Indeed, the community of the outer nearshore zone showed an increase in abundances and taxa numbers, and a decrease in diversity indices. Similar to an earlier study in the offshore zone of the nearby island of Norderney (Dörjes, 1976), the community was dominated by three characterizing species: Lanice conchilega, Spiophanes bombyx, and Magelona jonstonii. Information on the ecology of Magelona johnstonii and Spiophanes bombyx is limited. Studies showed that Magelona johnstonii and Spiophanes bombyx mainly occur in areas of fine sand. They are interface feeding polychaetes present in several subtidal areas of the North Sea (Fiege et al., 2000; Kröncke et al., 2011; Markert et al., 2015). All species probably benefit from food supply via nutrient discharge mediated growth of phytoplankton and microphytobenthos in the nearshore zone (Beck et al., 2017).

Lanice conchilega is a tube building polychaete colonizing both intertidal and subtidal sediments down to about 100 m depth (Dittmann, 1999; Ropert and Dauvin, 2000; Callaway, 2003). It occurs in high densities of several thousand individuals per m² on all coasts of Europe and in both the Atlantic and the Pacific (Holthe, 1977). Maximum abundances of up to 20,000 individuals per m² were reported from shallow subtidal areas of the German coast (Callaway, 2003). Tubes of Lanice conchilega protrudes 1 to 4 cm above the sediment surface and consist of smaller particles attached to an inner thin organic layer (Ziegelmeier, 1952; Dittmann, 1999). Studies on the ecology of Lanice conchilega showed that in both inter- and subtidal areas aggregations of tubes provide niches for a generally more species rich and abundant benthic faunal community (Rabaut et al., 2007; van Hoey et al., 2008). Lanice conchilega is, thus, often referred to as an important ecosystem engineer (Callaway, 2003; Rabaut et al., 2007; van Hoey et al., 2008; Callaway et al., 2010; de Smet et al., 2013). High abundances of Lanice conchilega increase the habitat quality in otherwise uniform habitats, and can increase sediment stability in exposed areas, resulting in a higher survival of the surrounding benthic species (Callaway, 2003; Rabaut et al., 2007; van Hoey et al., 2008; Callaway et al., 2010; de Smet et al., 2013). Indeed, the finding of a total number of 85 taxa, observed for the whole system, is surprisingly high compared to those known from literature. Here a total maximum of 25 taxa in a dissipative beach system were reported (McLachlan and Dorvlo, 2005).

The tubes of Lanice conchilega can also be detected by hydroacoustic measurements (Galloway and Collins, 1998; Ellingsen et al., 2002; Freitas et al., 2003; Wienberg and Bartholomä, 2005; Markert et al., 2013, 2015). Bartholomä (2006) stated that the occurrence of benthic faunal patterns covering a few m² immediately produces changes in backscatter signals. These changes correspond to shifts in the acoustic class assignment. For the outer nearshore zone of Spiekeroog we detected variation within the hydroacoustic signals by the backscatter measurements, becoming more obvious toward the offshore zone. In the offshore zone also current velocities of the BBL were lower compared to the nearshore zone. The community structure revealed again an increase in taxa number and in the total mean abundances. Eckman et al. (1981) showed that a decrease of the current velocities can be caused by an increase in habitat structure by tubes of Lanice conchilega; they cause deflections of the earlier laminar flow. The high abundances found in the offshore zone, thus, might have caused the decrease in the current velocities of the bottom layers.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.